The procedure works in the following way: Start with an initial guess for magnetic island widths for the resonant helical harmonics being considered. Determine the effect of each magnetic island on the background current density and pressure profile near each mode rational surface. The background axisymmetric current density and pressure profile and flux surface shapes (such as elongation and triangularity) can be taken, for example, from the evolving solution in an integrated modeling code. A solution is found for each harmonic of the perturbation from the set of coupled differential equations (7) and (8) together with ancillary equations (9) and (10) for and . The Hamada coordinate metric elements (the geometric terms in Eqs. (11-13)) are used to compute the harmonic mixing of the perturbed magnetic field components . An adaptive ordinary differential equation solver  can be used to integrate the equations for and from the magnetic axis to the plasma boundary or perfectly conducting wall. The island widths are adjusted iteratively until all the solutions match the boundary conditions at the wall. The width of each magnetic island acts like a nonlinear eigenvalue for the two point boundary value problem in which the perturbed magnetic field radial component vanishes at the magnetic axis and at a perfectly conducting wall at the plasma edge. For a self-consistent treatment, this algorithm is used in an integrated modeling code and the transport is enhanced across each island, as described in the next section.
At each stage in the iteration, the amplitudes of the perturbed magnetic field harmonics are adjusted to be consistent with the widths of the corresponding magnetic islands, using Eq. (28). In most other treatments of tearing modes, the differential equation for the perturbation is integrated from the magnetic axis outward to the island, and then from the wall inward to the island, and the solution is matched across the island. When there are multiple harmonics that are mixing with one another (through Eqs. (11-13)), it is simpler to integrate all the harmonics together from the magnetic axis to the outer boundary condition at the wall. With finite island widths and the resulting modifications of the profiles described in section 1.3 above, an adaptive ODE solver  can integrate all the harmonics from the magnetic axis through the islands to the outer boundary condition.
The algorithm described above applies to an arbitrary number of helical harmonics in any background axisymmetric geometry with arbitrary cross-section, aspect ratio, and plasma . The algorithm is simple enough to be used in an integrated predictive modeling code. Island widths change as the plasma profiles evolve. The algorithm does not apply to tearing modes that do not saturate in time, such as the tearing mode that drives sawtooth crashes.