Initial value methods. Boundary value problems were long considered some sort of a subclass of initial value problems (IVPs) whereby special care has to be taken of the initial conditions in order to get things right at the other end. This opinion has changed and IVPs are considered to be actually a special and in some sense relatively simple subclass of BVPs. The fundamental difference is that for IVPs one has the complete information about the solution at one point, i.e., the initial point, so some marching algorithm which is always local in nature may be considered. For BVPs, on the other hand, no complete information is available at any point, so the end points have to be connected by the solution algorithm in a global way. Only after stepping through the entire domain the solution at any point can be determined. In spite of the fundamental difference of the problem nature methods suited for initial values can also be applied to BVPs. The simplest initial value method is the shooting method. Values for all of the missing dependent variables at one boundary are chosen. These values must, of course, be consistent with any BCs for that boundary, but otherwise they are arranged to depend on arbitrary free parameters whose values are initially ``randomly'' chosen. Then the ODE system is integrated by any suitable initial value method, arriving at the other boundary. In general, there will be discrepancies from the desired boundary values there. However, by adjustment of the free parameters at the starting point the discrepancy can be zeroed. The BVP is thus transferred into a multidimensional root-finding problem. In case of linear ODEs a linear algebraic system has to be solved--a situation that applies to our problem. Unfortunately, this simple algorithm suffers from potential numerical instabilities. More advanced methods like reduced superposition, multiple shooting, reorthogonalization, and decoupling moderate these instabilities at the cost of increased numerical demands. The main advantage of initial value methods is the possibility of a very memory-saving implementation besides their conceptual simplicity since they somehow exploit the local nature of IVPs. In our case the number of unknowns is extremely high, the vertical variation of 4000 and more Fourier coefficients has to be determined. Hence, we implemented a very advanced IVP that follows the stabilized march algorithm derived in [200, pp. 155-164]. This algorithm has also successfully been implemented in the general purpose code MUSL [201] available through NETLIB.d Our implementation is tailored to the specific structure of the ODE (6.32) to enhance the performance; the core of the implementation, however, was taken from the MUSL code.
Finite-difference methods. Finite-difference or relaxation methods use a different approach. The ODEs are replaced by finite-difference equations on a mesh of points that covers the entire range of the integration interval. A trial solution of the values of the dependent variables at each mesh point is first guessed and thus generally does not satisfy the desired finite-difference equations, nor does necessarily satisfy the required BCs. Subsequent iterations also called relaxations try to adjust all the values on the mesh so as to bring them into successively closer agreement with the finite-difference equations and, simultaneously, with the BCs. Obviously the success of the relaxation depends on the quality of the guess, i.e., the closer it is to the true solution the better the performance. The number of values that have to be adjusted thus equals the number of unknowns multiplied by the number of mesh points. In our case this corresponds to roughly 4 105 values, i.e., 4000 unknowns and 100 mesh points. This number is far too big to be solved on today's workstations. Hence relaxation techniques are not suited for our problem. The dramatically high numerical demands result from the global solution strategy. In some cases--provided that the size of the problem is not too large--relaxation methods are the only means to solve the BVP since they work better than initial value methods when the BCs are especially subtle, or where they involve complicated algebraic relations that cannot easily be solved in closed form. In our case the BCs are simple linear relations and thus very easily evaluated. Furthermore, relaxation works better when the solution is smooth and not oscillatory. Such oscillations would require many grid points for accurate representation. Additionally the number and position of required points is not known a priori. Here lies another strength of initial value methods since they can employ variable step-size integrators to adjust naturally to the solution peculiarity. However, relaxation methods are advantageous in case of extraneous solutions which, while not appearing in the final solution satisfying all BCs, may wreak havoc on the initial value integrations required by shooting. A typical case is that of trying to maintain a dying exponential simultaneously with growing exponentials. The ODE system is then said to be stiff. Unfortunately, exactly this situation occurs in our problem since the waves traveling downwards the resist are damped whereas the waves traveling upwards the resist seem to be growing if one looks in the same direction both times, e.g., into the resist. Hence no simple shooting can be used in our case, and even the advanced techniques used can moderate such problems only up to a certain amount. The limitations of our algorithm thus caused are described in Section 6.5. Relaxation methods would overcome them but are prohibitively computation- and memory-intensive.