Solutions are obtained to large overdetermined systems of equations. Both nonlinear and linear systems are considered. The nonlinear system represents a dipole model of the earth's geomagnetic field, which is generated from spherical harmonic coefficients. This system of 64 unknowns and 1836 equations is solved by a maximum neighborhood method, which is an optimum interpolation between the well known Taylor's series and steepest descent methods. The original given values of the generated field are as large as 60,000 gamma, whereas a rms residual of 27.9 gamma is obtained with 173 iterations. The linear system of equations represents dipole changes required to account for the earth's secular change field which is generated from spherical harmonic coefficients. The dipole parameters computed from the nonlinear model are used as input parameters. The system contains 64 unknowns and 612 equations and is solved using a Chebyshev polynomial iterative method. These results are compared to results obtained by a direct solution of the normal equations of the system and results obtained by a pseudo-inverse method using a modified Gram-Schmidt factorization. Although the latter two methods give smaller rms values than the iterative method, the results of the iterative method are more reasonable in view of known properties of the results. The generated field has a rms value of 45 gamma per year. An rms residual of 2.5 gamma per year was obtained after 25,000 iterations.