We have developed a three-dimensional (3-D) algorithm for computing the controlled-source audio-frequency magnetotelluric response which is capable of handling local inhomogeneities in the electrical conductivity and magnetic permeability distribution. It is a difference equation algorithm that is based on the integral forms of Maxwell’s equations. To avoid a spatial singularity being introduced by the source, we used a secondary electric field approach. We adopt a staggered-grid finite difference method to disperse the secondary electric field equation which is simple to realize and can calculate relatively complex models. The final sparse equation system after dispersion is solved using an ILU preconditioning quasi-minimal residual (QMR) method. The code is validated for nonmagnetic and permeable conductive structures by two synthetic models. One of them is a two-dimensional (2-D) nonmagnetic conductive geoelectric model; we validate our 3-D computation result by comparing it to the result of 2-D CSEM finite element forward program. The other one is a 3-D permeable geoelectric model; we validate our 3-D computation result by using the COMSOL Multiphysics software. We also designed two theoretical 3-D geoelectrical models to analyze the effects of magnetic permeability heterogeneity. Numerical simulation results well reveal the influence of magnetic permeability heterogeneity on high-resistance and low-resistance models.