An optimal algorithm for certain boundary value problem

The O(h4) finite-difference scheme for the second derivative u″(x) leads to a coherent pentadiagonal matrix which is factorized into two tridiagonal matrices. This factorization is used to derive an optimal algorithm for solving a linear system of equations with the pentadiagonal matrix. As an application, a nonlinear system of ordinary differential equations is approximated by an O(h4) convergent finite-difference scheme. This scheme is solved by the implicit iterative method applying the algorithm at each iteration. A Mathematica module designed for the purpose of testing and using the method is attached.