INTERNATIONAL JOURNAL OF RESEARCH AND SCIENTIFIC INNOVATION (IJRSI)
ISSN No. 2321-2705 | DOI: 10.51244/IJRSI |Volume XII Issue X October 2025
Page 3873
www.rsisinternational.org
Block Optimized Hybrid Methods for Integrating Singular Second
Order Ordinary Differential Equations
Utalor I. Kate
1*
, Oyowei, E. Augustine
2
, OKafor M. Folakemi
1
, Ajie I. James
1
1
Department of Mathematics programme, National Mathematical Centre, Kwali Abuja.
2
Department of Statistics programme, National Mathematical Centre, Kwali Abuja.
*Corresponding Author
DOI: https://dx.doi.org/10.51244/IJRSI.2025.1210000334
Received: 10 November 2025; Accepted: 16 November 2025; Published: 22 November 2025
ABSTRACT
In this work, a block methods with characteristics of LMF are derived, analyzed and numerically applied to
solve singular Initial/Boundary value problems. It was done by applying shift operator to two linear multi-step
formula and combined with Optimize hybrid set of formula which are developed at the the first sub-interval to
circumvent the singularity at the left end of the integration interval.The mathematical derivation of the
proposed methods is based on method of undetermined coefficients where the coefficient in our Linear Multi-
step Formulas (LMF) are determined. The fundamental properties of the proposed scheme are analyzed.
Finally, the numerical implementation of the method are done on some singular I/B value problem which
demonstrate the accuracy and validity of the suggested technique when compared to various strategies
available in the current literature.
Keywords: One-block methods; shift operator, undetermined coefficients, Lane-Emden-type equation,
singular Initial/Boundary value problems (SIBVP).
INTRODUCTION
The goal of this research paper was to find a reliable numerical approach for the solution of the singular
Initial/Boundary value problem (SIBVP) of Lane–Emden equations of the form:
′′
′ (1)
subject to the boundary conditions
′
′
(2)
where
are real constants, is continuous real function. The Existence and
uniqueness of the solution to the problem (1) subject to any boundary conditions have been rigorously
determined by Zhang [24].
Second-order singular boundary value problems are commonly encountered in several areas of applied
mathematics, physics and engineering, such as mathematical modeling, chemical kinetics, astrophysics,
catalytic diffusion reactions [25] and among others. Researchers in various fields such as applied sciences and
engineering have shown significant interest in solving equations (1) by trying to find better and efficient
methods.The problem under consideration becomes one of the most complex problems to solved analytically,
due to the nonlinear properties of (1) and the singularity arising at the point t = 0, called singular point.In order
to overcome these challenges and obtain meaningful solutions, numerical methods have emerged as crucial
INTERNATIONAL JOURNAL OF RESEARCH AND SCIENTIFIC INNOVATION (IJRSI)
ISSN No. 2321-2705 | DOI: 10.51244/IJRSI |Volume XII Issue X October 2025
Page 3874
www.rsisinternational.org
tools, example of such methods include the finite difference methods proposed in Kumar [13] and Pandey
[16], the spline methods discussed in Caglar et al.,[8], Kadalbajoo et al. [12], the approximation methods
introduced in (Allouche et al. [5], and Aydinlik et al., [6], the high-order compact finite differences method in
Malele et al.,[15] and among others.
In recent time, researchers has employ the use of the optimization technique by Ramos et al.,[17], Asifa et al.,
[1] and Rajat et al. [19] and non-optimization by Jator,[11] and Anake et al., [4] in solving general second
order problem. .
The focus of this paper is block methods which posses good stability properties for solving differential
equations. They are constructed using two different LMF with aid of shift operator, which are combined with
optimized hybrid set of formula called ad-hoc method that is applied only to the first sub interval due to the
singularity at t = 0. we aimed to obtain the optimization formulation ad-hoc developed in Utalor et al., [23] to
further improve and check their performance. In this way, we obtain a scheme capable of solving the problem
posed effectively.
The present work is outlined as follows. In Section 2, we present the KSPHT method for solving SBVPs. The
characteristics of the developed formulas are analyzed in Section 3. In Section 4,shows the Implementation of
the method.We present the numerical results of some Test problems to show the efficiency and reliability of
the proposed technique in Section 5. Conclusions are outlined in Section 6.
CONSTRUCTION OF THE METHOD
We approximate the exact solution of (1) in the partition
of the
integration interval
, with constant step size
by a self-starting
block method. The continuous coefficients (
) of the composing LMF are
determined by imposing order condition on linear muti-step formula (L.M.F) and using the method of
undetermined coefficients developed in Ajie et al.,[2, 3] and Brugnano et al., [ 7 ]. See [23] for how the the
self-starting block methods (KSPHT) for the main method are obtained, and also the 2 off-steps of non-
optimization formulas to circumvent the singularity.
First, reformulate the equation (1)
″
′
ℎ
′
′
(3)
Thus, the singularity is transferred to the function .This block method cannot be used directly for solving a
BVP problem in (1) because it is not possible to evaluate
, since there is a singularity at
.To overcome this drawback, we develop a set of multi-step formulas to be applied at the first sub-interval
′
ℎ
(4)
with the purpose of specifically avoiding the use of
. as a result of it, the method will have main
formula and also the Formulas to Circumvent the Singularity.
MAIN FORMULAS
Let us consider the Linear Multi-step method(LMM) of the form
′
ℎ
ℎ
INTERNATIONAL JOURNAL OF RESEARCH AND SCIENTIFIC INNOVATION (IJRSI)
ISSN No. 2321-2705 | DOI: 10.51244/IJRSI |Volume XII Issue X October 2025
Page 3875
www.rsisinternational.org
′
ℎ
ℎ
using definition of order, This leads to the following matrix equation:
′
For
In equation (6) when is solved by Mathematica software package method to obtain the value of the
continuous coefficient
′
ℎ
′
ℎ
ℎ
and its derivative as
′
′
ℎ
′
ℎ
ℎ
Evaluating (7) at the points t=3 gives the method and its derivative.
Applying the theory in Utalor et al., [23] on the method and its derivative , the coefficients of the resultant
block method after the shift operator application in vector form are below
INTERNATIONAL JOURNAL OF RESEARCH AND SCIENTIFIC INNOVATION (IJRSI)
ISSN No. 2321-2705 | DOI: 10.51244/IJRSI |Volume XII Issue X October 2025
Page 3876
www.rsisinternational.org
-
-
-
-
(8)
OPTIMIZATION FORMULAS TO CIRCUMVENT THE SINGULARITY (ONE-STEP
METHOD WITH TWO OPTIMIZE POINTS
Considering different intermediate points using undetermined coefficient .This off-step point are gotten by
minimization of Local Truncation Error of the intermediate points of the main formula at the grid points, as to
circumvent the singularity at the left end of the integration interval, as a result we will developed a set of multi-
step formulas specially designed for the sub-interval
t
0
,t
1
, where the value f
0
is absent,
Let us consider the Hybrid Linear Multi-step method(HLMM) of the form
′
ℎ
This leads to the following matrix
equation:
′
Applying equation (9), two two different intermediate points are introduced,
, Where
Equation (10) is solved by Mathematica
software package method to obtain the value of the unknown parameters
expressed as functions of t (whose expressions are not included), and can be written as
ℎ
′
ℎ
′
ℎ
Evaluating (11) at the points gives the continuous form of the method, which implied that
ℎ
′
ℎ
INTERNATIONAL JOURNAL OF RESEARCH AND SCIENTIFIC INNOVATION (IJRSI)
ISSN No. 2321-2705 | DOI: 10.51244/IJRSI |Volume XII Issue X October 2025
Page 3877
www.rsisinternational.org
ℎ
′
ℎ
ℎ
′
ℎ
The first derivative of equation (9) with respect to t gives
ℎ
′
′
ℎ
Evaluating (13) at the points gives the addition method, which implied that
′
′
ℎ
′
′
ℎ
′
′
ℎ
In order to determine appropriate values for r,s, we choose to optimize the local truncation errors in the main
formulae (13 and 14) respectively. which is obtained after expanding in Taylor series around t
n
, which results
in
)(0
72
))())62(21(
));((
)(0
360
))())41(552(
));((
6
5)5(
1
6
5)5(
1
h
htysrs
htyhL
h
htysrs
htyL
n
n
n
n
(15)
Equating the principal term of this error to zero in each term in (15) , that is the coefficients of h
5
in the above
formulae, we obtain the system
and solving (16) for r and s ,we get the value as
and thus, there is a unique solution with the constraints 0<r<s<1 .
INTERNATIONAL JOURNAL OF RESEARCH AND SCIENTIFIC INNOVATION (IJRSI)
ISSN No. 2321-2705 | DOI: 10.51244/IJRSI |Volume XII Issue X October 2025
Page 3878
www.rsisinternational.org
Note that we have six unknowns (
′
), to get a one-step hybrid block method for solving the
SBVP problem, we need to complete the above formulas. For that we consider the evaluation at y
r
, y
s
and its
derivative. Considering the values in the block method results to be the following system of six equations
ℎ
′
ℎ
ℎ
′
ℎ
ℎ
′
ℎ
′
′
ℎ
′
′
ℎ
′
′
ℎ
ANALYSIS OF THE METHODS
Order and error Constants of the Methods [14]
The linear difference operator associated with the block (8) is defined
ℎ
ℎ
ℎ
ℎ
Expanding (8) using Taylor series ,we obtained
ℎ
ℎ
′
ℎ
″
ℎ
ℎ
So that
ℎ
ℎ
ℎ
ℎ
ℎ
Here p is the order and
is the error constant (Lambert, 1973). The following table shows the error constant
Table 1. Error constants of the composing LMFs
formulae
formulae
where the formulae , represent
, it follows that for all the formulae, the order
Since the order of the formulas is greater than one, they are consistent. For the ad- hoc formulas used for the
first step, it is easy to see that they are also consistent
Zero Stability
Definition 1:The implicit block method (8) is said to be zero stable if the roots z
s
, s=1,…,n of the first
characteristic polynomial ρ
z
, defined by
L
INTERNATIONAL JOURNAL OF RESEARCH AND SCIENTIFIC INNOVATION (IJRSI)
ISSN No. 2321-2705 | DOI: 10.51244/IJRSI |Volume XII Issue X October 2025
Page 3879
www.rsisinternational.org
satisfies /z
s
/≤1 and every root with /z
s
/=1 has multiplicity not exceeding two in the limit as h→0 [14, 25].
Using the definitions, the method in (8) may be rewritten in a more appropriate vector form to study zero-
stability as
ℎ
ℎ
The same procedure are done for the ad-hoc formulas used for the first step,(whose expressions are not
included), its was proofed to be zero stable and have higher order more than the the non- optimization formula
due to the optimize strategy done.
Convergence
Theorem (1): Consistency and zero stability are sufficient condition for linear multistep method to be
convergent [14]. Since the method (8) is consistent and zero stable, it implies the method is convergent for all
point .
Implementation
The derive KSPHT which include the ad-co formula are combined and applied as block form. The solutions are
considered in the interval
where N is the number of blocks. The formulas are
written as G(y) = 0 with the following unknown values to be obtained and the unknowns are expressed as
′
′
′
The resulting system is solved using Newton-Raphson iteration given as
change_in_y =
;
y_new = change_in_y+(y_old);
y_old = y_new;
where J is the Jacobian matrix of G. The following Taylor’s approximations are considered as starting values
ℎ
′
ℎ
′
′
ℎ
NUMERICAL ILLUSTRATIONS
Numerical examples are presented in this section to show the efficiency of the developed methods, K-step pair
of hybrid techniques (KSPHT) which include 1S2OP(one step, Two optimized points). The accuracy is
measured by using the following formulas:
, where denotes the absolute error at the
INTERNATIONAL JOURNAL OF RESEARCH AND SCIENTIFIC INNOVATION (IJRSI)
ISSN No. 2321-2705 | DOI: 10.51244/IJRSI |Volume XII Issue X October 2025
Page 3880
www.rsisinternational.org
considered node,
and
are exact and approximate solutions of the problems, respectively.
The following notations are used in the tables when presenting the results:
Block Hybrid Methods which include non-optimazition of (KSPHT) 1S2HP , [23]. TWS-Taylor wavelet
solution [10], AADM - Advanced Adomian decomposition method, [22], MLMF - Modified Linear Multistep
Formulas [18]. The computed results for the three problems using the methods proposed are presented in tables
and graphically.
Problem 1.
Consider the following physical model SBVP problem of the isothermal gas sphere equilibrium, as described
in Gumgum, and Umesh et al., [10, 22]:
′′
′
′
(4.1.2)
The equation arise in the study of stellar structure where m=5. It exact solution is
. The
example is solved within the interval [0, 1] over ten (10) iterations and the results are compared both with the
other numerical results and the exact solution to show the efficiency and validity of the method.
Table 2. Comparison of absolute errors of Problem 1 obtained using KSPHT (1S2OP)
x
Exact
1S2OP
1S2HP
AADM
TWS
0.1
0.998337488459583
2.3024615e-8
7.78932376e-7
1.65000e-6
6.46000 e- 6
0.2
0.993399267798783
2.2271241e-8
7.5386121e-7
6.63000e-6
6.30000e- 6
0.3
0.985329278164293
2.0803594e-8
7.07239719e-7
1.59000e-6
6.05000e- 6
0.4
0.974354703692446
1.8716236e-8
6.3830484e-7
1.53000e-6
5.70000 e- 6
0.5
0.960768922830523
1.6096693e-8
5.49157323e-7
1.44000-6
5.30000 e- 6
0.6
0.944911182523068
1.3074469e-8
4.44458267e-7
1.34000e-6
4.84000 e- 6
0.7
0.927145540823120
9.805192e-9
3.30551613e-7
1.10000e-6
4.33000 e- 6
0.8
0.907841299003204
6.44704e-9
2.14367027e-7
9.58000e-7
3.86000e- 6
0.9
0.887356509416114
3.141175e-9
1.02404574e-7
7.30000e-7
3.24000 e- 6
1
0.866025403784439
0
0
1.8900e-14
1.45000e- 13
INTERNATIONAL JOURNAL OF RESEARCH AND SCIENTIFIC INNOVATION (IJRSI)
ISSN No. 2321-2705 | DOI: 10.51244/IJRSI |Volume XII Issue X October 2025
Page 3881
www.rsisinternational.org
Figure 1. Plots of exact and KSPHT solution for Problem (1). It show good agreement between the numerical
and exact solutions.
PROBLEM 2
Consider the nonlinear heat conduction model of the human head, ′′
′
′
′
Where . The above nonlinear SBVP is discussed by Duggan R and Goodman A [19] as a heat
conduction model in the human head. However, The general analytical solution of problem is unknown [22].
Table 3. Comparison of KSPHT and the exact solution on Problem 2
x
1S2OP
1S2HP
MLMF
AADM
0
0.270029664529952
0.270029706093745
0.2700296478967
0.2700296466
0.1
0.268756917547862
0.268756958875822
0.2687569006296
0.2687568993
0.2
0.264932833272991
0.264932875883028
0.2649328175383
0.2649328162
0.3
0.258539803870774
0.258539847740317
0.2585397893815
0.2585397881
0.4
0.249548193506706
0.249548238415984
-
0.2495481789
0.5
0.237915899626578
0.237915945201108
-
0.2379158863
0.6
0.223587718025932
0.223587763723984
-
0.2235877058
0.7
0.206494492698337
0.206494537783022
0.2064944830238
0.2064944817
0.8
0.186552022693196
0.186552066192176
0.1865520141667
0.1865520128
0.9
0.163659688980158
0.163659729631207
0.1636596815804
0.1636596802
1
0.137698752907342
0.137698789086042
0.1376987466136
0.1376987453
INTERNATIONAL JOURNAL OF RESEARCH AND SCIENTIFIC INNOVATION (IJRSI)
ISSN No. 2321-2705 | DOI: 10.51244/IJRSI |Volume XII Issue X October 2025
Page 3882
www.rsisinternational.org
Figure 2. Plots of exact and KSPHT solution for Problem 2.
Problem 3.
The following model which corresponds to the reaction–diffusion process in a spherical permeable catalyst as
reported in Allouche [5], Utalor et al. [23]
′′
′
′
its analytical solution for n = 1, is given by
ℎ
ℎ
ℎ
whererepresents the Thiele modulus. The value ofis determined by
.
Table 4. Comparison of absolute errors of Problem 3 using KSPHT (1S2OP) within [0,1] over ten (10)
iterations
x
Exact
1S2OP
1S2HP
0.1
0.070225439227791
1.16507398e-7
3.36833567e-7
0.2
0.079188028691280
1.02523148e-7
3.67734455e-7
0.3
0.095650823305130
9.3155702 e-8
4.19915345e-7
0.4
0.122193513582708
8.5924964 e-8
4.92803425e-7
0.5
0.163071231929978
7.8751406e-8
5.83054692e-7
0.6
0.225009916448920
6.9694737e-8
6.80232762e-7
0.7
0.318481161564393
5.7022468e-8
7.58832348e-7
0.8
0.459715910279232
3.9708057e-8
7.63576405e-7
0.9
0.673870380204315
1.8740297e-8
5.82698005e-7
INTERNATIONAL JOURNAL OF RESEARCH AND SCIENTIFIC INNOVATION (IJRSI)
ISSN No. 2321-2705 | DOI: 10.51244/IJRSI |Volume XII Issue X October 2025
Page 3883
www.rsisinternational.org
1
1
0
0
To analyze the impact of the Thiele modulus () on the concentration profile (y(x)),
we also considered other values of and . Figure 3 displays the numerical outcomes for various values of
and . We observed that in Figure 3, the concentration profile increases when diminishes.
Figure 3. Plots of exact and KSPHT solution for Problem 3.
DISCUSSION OF RESULTS
The results obtained from the three test problems are summarized in Tables 2–4 and Figures 1–3. In Table 2,
the solution of Problem 1 obtained using the proposed methods (KSPHT), which include 1S2OP, at the points
x = 0(0.1)1.0, is compared with the results of Utalor et al. (2025), Umesh et al. (2021), and Gumgum (2020).
Overall, the KSPHT method based on the optimization technique demonstrates superior performance compared
to other existing methods, as shown in column three of Table 2. Table 3 and Figure 2 present the comparison
of approximate solutions for Problem 2. The proposed results show very good agreement—up to seven to eight
decimal places—with those obtained using the Modified Linear Multistep Formulas (Olabode et al., 2024), the
Advanced Adomian Decomposition Method (Umesh et al., 2021), and the Block Hybrid Methods (Utalor et al.,
2025). For Problem 3, the results displayed in Table 4 and Figure 3 indicate that the proposed methods
(KSPHT), particularly the 1S2OP variant, outperform the Block Hybrid Methods of Utalor et al. (2025).
Figures 1 and 3 further show that the numerical solutions are in close agreement with the exact solutions. In
general, the proposed methods compare favorably with existing approaches in the literature, demonstrating
improved accuracy and efficiency despite differences in formulation.
CONCLUSION
A new approach for constructing self-starting block methods for solving second-order singular boundary value
problems (SIBVPs) has been presented. The strategy involves applying a shift operator to two distinct linear
multistep formulas and combined with an optimized hybrid set of formulas developed over the first sub-
interval. The continuous coefficients of the linear multistep methods were obtained using the method of
INTERNATIONAL JOURNAL OF RESEARCH AND SCIENTIFIC INNOVATION (IJRSI)
ISSN No. 2321-2705 | DOI: 10.51244/IJRSI |Volume XII Issue X October 2025
Page 3884
www.rsisinternational.org
undetermined coefficients. To demonstrate the effectiveness of the proposed methods, three real-world model
problems from the literature were solved. The numerical results show that the proposed methods produce
highly accurate solutions with smaller errors when compared to existing techniques. Furthermore, the results
indicate that the optimized selection of off-step points yields superior performance compared to non-optimized
approaches, confirming the efficiency and robustness of the developed methods.
REFERENCES
1. Asifa, T., Sania, Q., Amanullah , S., Evren , H., & Asif, A. (2022). A new continuous hybrid block
method with one optimal intra step point through interpolation and collocation Fixed Point . Theory
Algorithms Sci Eng , 22.
2. Ajie, I. J., Durojaye, M.O, Utalor, I. K., & Onumanyi, P. A. (2020). Construction of a family of stable
one- Block methods Using Linear Multi- Step Quadruple. Journal of Mathjematics., e-ISSN: 2278-
5728
3. Ajie, I. J., Utalor, I. K., & Onumanyi, P. A. (2019). A family of high order one-block methods for the
solution of stiff initial value problems. Journal of Advances in Mathjematics and Computer Science.,
31(6), 1–14.
4. Anake , T. A., Awoyemi, D. O., & Adesanya , A. O. (2012). One-step implicit hybrid block method
for the direct solution of general second order ordinary differential equations,. Int. J. App. Math. ,
42(4), 224–228.
5. Allouche, H.; Tazdayte, A. (2017). Numerical solution of singular boundary value problems with
logarithmic singularities by Padè approximation and collocation methods. J. Comput. Appl. Math. ,
311(324–341.).
6. Aydinlik, S.; Kiris, A. (2018). A high-order numerical method for solving nonlinear Lane-Emden-type
equations arising in astrophysics. Astrophys. Space Sci., 363, 264.
7. Brugnano, L., Trigiante, D. (1998). Solving Differential Problems by Multistep Initial and Boundary
Value Methods. Amsterdam: Gordon and Breach Science Publishers.
8. Caglar, H.; Caglar, N.; Ozer, M. (2009). B-spline solution of non-linear singular boundary value
problems arising in physiology. Chaos Solitons Fract., 39(1232–1237.).
9. Duggan R & Goodman A. (2017). Pointwise bounds for a nonlinear heat conduction model of the
human head. Bulletin of Mathematical Biology,Springer Science and Business Media LLC, 48(2),
229-236. DOI:10.1016/s0092-8240(86)80009-x
10. Gumgum, S. (2020). Taylor wavelet solution of linear and non-linear Lane-Emden equations. Appl.
Numer. Math,158, 44–53.
11. Jator, S.N.: On a class of hybrid methods for y = f (x,y,y ). Int. J. Pure Appl. Math. 59(4), 381– 395
(2010)
12. Kadalbajoo, M.K.; Kumar, V. (2007). B-Spline method for a class of singular two-point boundary
value problems using optimal grid. Appl. Maths. Comput. , 188(1856–1869.).
13. Kumar, M. A . (2003). difference scheme based on non-uniform mesh for singular two-point boundary
value problems. Appl. Math. Comput. , 136(281–288.).
14. Lambert J. D.(1973), Computational methods in Ordinary Differential Equations, John Wiley and
sons, New York.
15. Malele, J.; Dlamini, P.; Simelane, S. (2023). Solving Lane–Emden equations with boundary
conditions of various types using high-order compact finite differences. Appl. Math. Sci. Eng. , 31,
2214303.
16. Pandey, R.K. (1997). A finite difference methods for a class of singular two point boundary value
problems arising in physiology. Int. J. Comput. Math., 65(131–140.).
17. Ramos, Z. Kalogiratou, Th. Monovasilis and T. E. Simos, An optimized two-step hybrid block method
for solving general second order initial-value problems, Numerical Algorithm, Springer 2015. DOI
10.1007/s11075-015- 0081-8.
18. Olabode , B. T., Kayode , S. J., Odeniyan-Fakuade , F. H., & Momoh , A. L. (2024). Numerical
Solution to Singular Boundary Value Problems (SBVPS) using Modified Linear Multistep Formulas
(LMF). International Journal of Mathematical Sciences and Optimization: Theory and Applications,
INTERNATIONAL JOURNAL OF RESEARCH AND SCIENTIFIC INNOVATION (IJRSI)
ISSN No. 2321-2705 | DOI: 10.51244/IJRSI |Volume XII Issue X October 2025
Page 3885
www.rsisinternational.org
10(1), 34–52.
19. Rajat, S., Gurjinder, S., Higinio , R., & Kanwar , V. (2023). An efficient optimized adaptive step-size
hybrid block method for integrating w′′ = f (t,w,w′) directly. Journal of Computational and Applied
Mathematics, 114838, 420.
20. Rufai, M. A., and Ramos, H. (2021). Numerical Solution for Singular Boundary Value Problems
Using a Pair of Hybrid Nyström Techniques.Axioms 2021, 10, 202. Axioms , 10, 202.
21. Thula, K.; Roul, P. A High-Order B-Spline Collocation Method for Solving Nonlinear Singular
Boundary Value Problems Arising in Engineering and Applied Science. Mediterr. J. Math. 2018, 15,
176.
22. Umesh M, & Kumar, M. (2021). Numerical solution of singular boundary value problems using
advanced Adomian decomposition method. Eng. Comput., 37, 2853–2863.
23. Utalor, I. K., Ajie, I. J., Durojaye, M.O, (2025). Construction of Block Hybrid Methods for the
Solution of Singular Second Order Ordinary Differential Equation. Asian Research Journal of
Mathematics,21(5),23–37. https://doi.org/10.9734/arjom/2025/v21i5923
24. Zhang, Y. (1993). Existence of solutions of a kind of singular boundary value problem. Nonlinear
Anal. Theory Methods Appl. , 21(153–159.).
25. Henrici, P.(1962). Discrete Variable Methods in ODE. New York: John Wiley & Sons.
