We have investigated different approaches to solve magnetotelluric forward modeling problems. Besides classical direct electromagnetic (EM) formulation, ungauged and gauged (Lorenz, Coulomb, and axial) vector and scalar potential formulations are discretized using the finite-difference (FD) numerical solution technique. Linear matrix equations obtained from FD discretization for each EM field formulation are solved by using the conjugate-gradient iterative solution method. We compared FD solutions of each approach with respect to accuracy and speed. We found that all approaches generate accurate results. However, an ungauged approach gave solutions with a lower conjugate-gradient iteration and accordingly less computer time. All stiffness matrices arising from each formulation are examined in terms of their size and type. FD solution of the Coulomb-gauge and ungauged approaches generates larger stiffness matrices than the other approaches. Direct EM, ungauged, and axial-gauge approaches generate a symmetric stiffness matrix. We recommend the use of the ungauged approach in inversion algorithms, which gives faster forward solution, amongst other things.