The magnetic flux evolution problem in toroidal plasmas is formulated in a framework suitable for integrating externally imposed magnetic field components with internal components from bootstrap current and auxiliary current drive. The formulation is applicable to 3-dimensional (3-D) stellarator equilibria, and reduces to 2-D form for axisymmetric plasmas. Here the numerical implementation of this framework is described. Conservative integrations schemes, resolution close to the magnetic axis, and efficient methods for flux surface averaging are discussed. Results from the test code THRIFT (THRee dimensional Inductive Flux evolution in Toroidal plasmas) are used to illustrate numerical convergence properties for a low aspect ratio stellarator and the axisymmetric NSTX spherical torus.