High order fitted operator numerical method for self-adjoint singular perturbation problems (Q814747)

A numerical scheme is proposed for self-adjoint singular perturbation problems in conservation form \[ -\varepsilon (a(x)y')' + b(x)y = f(x), \quad x \in [0,1] \] subject to nonhomogeneous Dirichlet boundary conditions, for sufficiently smooth data satisfying \(a(x) \geq a > 0, b(x) \geq b >0\), and \(0 < \varepsilon \ll 1\). The same problem has been numerically solved by using a cubic spline difference scheme with exponentially fitting factor \textit{M. K. Kadalbajoo} and \textit{K. C. Patidar} [Int. J. Math. Math. Sci. 2003, No. 61, 3873--3891 (2003; Zbl 1045.65064)]. The same trick is adapted here just by replacing the cubic spline scheme by the Numerov scheme with the same fitting factor. The author is claiming that the present scheme is of fourth-order accurate for moderate values of the singular perturbation parameter \(\varepsilon\), and for very smaller values of \(\varepsilon\) the method is \(\varepsilon\)-uniform convergent with order two. But this is not clear from the theorem on the error estimate. Indeed, the standard Numerov's scheme (i.e., without any fitting factor) produces good approximate results on a uniform mesh for moderate values of \(\varepsilon\). From the numerical tables on the convergence rate it is very well clear that the method is not fourth-order convergent. Also, one can notice that for certain values of \(\varepsilon\) the rate convergence will be less than one, but there is no explanation about these strange behavior.