We present a MATLAB solver designed with special focus on the following class of nonlinear singular boundary value problems with a singularity of the first kind:

z'(t)=M(t)z(t)/(t-a)+f(t,z(t)), t∈ (a,b], B_az(a)+B_bz(b)=β,

However, this code can also be applied to treat the more general problem class,

z'(t)=F(t,z(t)), t∈ [a,b], B(z(a),z(b))=0,

(2)

The solver is based on collocation at either equidistant or Gaussian collocation points. For singular problems, the choice of equidistant nodes does not imply a loss of efficiency since no convergence order higher than the stage order can be expected in general. An asymptotically correct error estimate for the global error of the approximate solution is also provided. This estimate is a modification of an idea originally proposed by Zadunaisky and based on the defect correction principle. Particular attention was paid to the efficiency of the error estimation procedure. On the basis of this estimate, we succeeded to realize a mesh selection procedure which is robust with respect to the singularity. This means that if the solution is smooth close to the singular point the mesh remains coarse in this region, independently of the unsmoothness of the direction field. We discuss the details of the implementation in MATLAB 5.3 (with a preview of a MATLAB 6.0 version), describe the measures taken to optimize the performance of the code and demonstrate how the mesh adapting works.

References

[1] W. Auzinger, O. Koch, and E. Weinmüller. Efficient collocation schemes for singular boundary value problems. Submitted to Numer. Algorithms.

[2] R. Frank and C. Überhuber. Iterated Defect Correction for differential equations, Part I: Theoretical results. Computing, 20:207-228, 1978.

[3] P.E. Zadunaisky. On the estimation of errors propagated in the numerical integration of ODEs. Numer. Math., 27:21-39, 1976.