It is commonly believed that the divergence equations in the Maxwell equations are "redundant" for transient and timeharmonic problems, therefore most of the numerical methods in computational electromagnetics solve only two first-order curl equations or the second-order curl-curl equations. This misconception is the true origin of spurious modes and inaccurate solutions in computational electromagnetics. By studying the div-curl system this paper clarifies that the first-order Maxwell equations are not "overdetermined," and the divergence equations must always be included to maintain the ellipticity of the system in the space domain, to guarantee the uniqueness of the solution and the accuracy of the numerical methods, and to eliminate the infinitely degenerate eigenvalue. This paper shows that the common derivation and usage of the second-order curl-curl equations are incorrect and that the solution of Helmholtz equations needs the divergence condition to be enforced on an associated part of the boundary. The div-curl method and the least-squares method introduced in this paper provide rigorous derivation of the equivalent second-order Maxwell equations and their boundary conditions. The node-based least-squares finite element method (LSFEM) is recommended for solving the first-order full Maxwell equations directly. Examples of the numerical solutions by LSFEM are given to demonstrate that the LSFEM is free of spurious solutions.