We will only outline the principle of multiple shooting employing
reorthogonalization and decoupling. Additionally, compactification is
applieda technique that leads to formulae that can be evaluated recursively.
The storage demands of the algorithm are thus not increased in comparison
to the previous ones and the performance with respect to stability is not
degraded due to the separated BCs.^{g}
Both properties are of utmost importance for our application.
The description of this socalled stabilized march algorithm is divided
into three
different parts since the first and last subinterval have to be treated
individually.
A graphical illustration of the algorithm, especially the integration direction
and the multiple shooting technique, is shown in Figure 6.4.

First interval (p = 1). In the first interval z [z_{1}, h] we carry out Step 1 and Step 2 of the shooting algorithm using reduced superposition (cf. Table 6.3), i.e., we perform a QRfactorization of the boundary matrix to determine the initial values (h) = at the first starting point z_{0} = h, and integrated then the IVP
(6.46) 
Interior intervals (p = 2,..., N_{p}  1). For the interior intervals z [z_{p}, z_{p  1}] we first determine the initial values (z_{p  1}) at the starting point z_{p  1}. Simply using the integrated solution (z_{p  1}) from the previous interval would be identical to single shooting. A common choice is taking any orthogonal matrix independently from the previous interval. The algorithm is then referred to standard multiple shooting [200, pp. 146148] and has the big advantage that the integrations within the N_{p} subintervals can be performed parallel. However, all coefficient matrices have to be stored individually. After integration they are determined by matching the solutions at the shooting points z_{p}. To obtain a stable algorithm this has to be done simultaneously [200, pp. 149153], i.e., a linear algebraic equation system with N_{p} x N/2 unknowns and N_{s} righthand sides has to be solved. Although this equation system is bandlimited with bandwidth N/2 the storage demands are dramatically increased by the number N_{p} of subintervals. A dangerously appealing idea is to evaluate the relations recursively. The memory consumption would be the same as for single shooting with reduced superposition but the stability properties are almost lost. This method, usually called multiple shooting with compactification [200, pp. 153155], is therefore not suited for our application.
However, if the initial values are not independently chosen from the previous interval the algorithm thus obtained becomes stable and is said to employ marching techniques [200, pp. 155164]. In order to achieve stability, it is necessary to monitor the growth of the fundamental modes. Otherwise the growing modes would quickly get dominant and the linear independence would get lost numerically due to the finite numeric precision of computers. Such a monitoring is realized by reorthogonalization or decoupling that is done by QRfactorization of the integrated solution matrix (z_{p  1}) of the previous interval p  1, i.e.,
(6.47) 
(6.48) 
Last interval p = N_{p}. For the last interval z [0, z_{Np  1}] the QRfactorization of the integrated solution (0) is skipped. The coefficient matrix is calculated by enforcing the BCs valid at z = 0, i.e.,
(6.49) 
(6.50) 