Using a three-component magnetic field data set at over 100,000 satellite points previously compiled for spherical harmonic analysis, we have produced a continuously varying magnetization model for Mars. The magnetized layer was assumed to be 40 km thick, an average value based on previous studies of the topography and gravity field. The severe nonuniqueness in magnetization modeling is addressed by seeking the model with minimum root-mean-square (RMS) magnetization for a given fit to the data, with the trade-off between RMS magnetization and fit controlled by a damping parameter. Our preferred model has magnetization amplitudes up to 20 A/m. It is expressed as a linear combination of the Green's functions relating each observation to magnetization at the point of interest within the crust, leading to a linear system of equations of dimension the number of data points. Although this is impractically large for direct solution, most of the matrix elements relating data to model parameters are negligibly small. We therefore apply methods applicable to sparse systems, allowing us to preserve the resolution of the original data set. Thus we produce more detailed models than any previously published, although they share many similarities. We find that tectonism in the Valles Marineris region has a magnetic signature, and we show that volcanism south of the dichotomy boundary has both a magnetic and gravity signature. The method can also be used to downward continue magnetic data, and a comparison with other leveling techniques at Mars' surface is favorable.