Numerical Study of Some Nonlinear Boundary Value Problems (BVPs) with an Efficient Numerical Method

Md. Amirul Islam^1*, Nazmun Nahar², S. M. Kamal Hossain²

¹Department of Mathematics, Mawlana Bhashani Science and Technology University, Tangail, Bangladesh.
²Department of Mathematics, Uttara University, Dhaka, Bangladesh.

DOI: 10.4236/ajcm.2025.153015                

Abstract

This study presents a comprehensive numerical investigation of fourth-order nonlinear boundary value problems (BVPs) using an efficient and accurate computational approach. The present work focuses on a class of nonlinear boundary value problems that commonly emerge in scientific and engineering applications, where the underlying models are often governed by complex nonlinear differential equations. Due to the difficulty of obtaining exact analytical solutions for such problems, numerical techniques become essential for reliable approximation. In this work, the Finite Difference Method (FDM) is adopted as the core numerical tool due to its robustness and effectiveness in solving such problems. A carefully designed finite difference scheme is developed to discretize the governing fourth-order nonlinear differential equations, converting them into a system of nonlinear algebraic equations. These systems are subsequently solved numerically using Maple software as the computational tool. The article includes two illustrative examples of nonlinear BVPs to demonstrate the applicability and performance of the proposed method. Numerical results, including graphical representations, are provided for various step sizes. Both absolute and relative errors are calculated to assess the accuracy of the solutions. The numerical findings are further validated by comparing them with known analytical or previously published approximate results. The outcomes confirm that the finite difference approach yields highly accurate and reliable solutions.

Keywords

Nonlinear Boundary Value Problems (BVPs), Finite Difference Method (FDM), Numerical Solutions, Absolute Error, Relative Error, Numerical Examples

1. Introduction

Nonlinear boundary value problems (BVPs) involving ordinary differential equations (ODEs) play a vital role in both theoretical and applied mathematics. These problems frequently arise in fields such as physics, engineering, and computational science, where systems are modeled by nonlinear differential equations with constraints imposed at the boundaries of a domain. Examples include the bending of elastic beams, steady-state heat distribution, nonlinear electrical circuits, chemical reaction kinetics, and mechanical systems under external forces. Unlike initial value problems (IVPs), BVPs often require boundary conditions at multiple points, making traditional step-by-step methods like Euler or Runge-Kutta less suitable. Moreover, the nonlinear nature of these problems introduces additional challenges, such as the lack of closed-form analytical solutions, structural complexity, and sensitivity to initial guesses when iterative methods are used. These complexities highlight the necessity of accurate and efficient numerical techniques. Despite significant progress in numerical analysis, solving nonlinear BVPs remains a challenging task due to potential issues with stability, convergence, and computational efficiency. Among the many available methods, the Finite Difference Method (FDM) stands out due to its simplicity, flexibility, and effectiveness in discretizing differential equations. However, when applied to nonlinear problems, FDM must be carefully tailored to maintain accuracy and ensure numerical stability.

This study aims to investigate the numerical behavior and performance of FDM in solving selected nonlinear BVPs. By constructing a suitable finite difference scheme and applying iterative techniques, the study evaluates how effectively FDM can approximate solutions to the nonlinear boundary value encountered in applied mathematics and engineering.

The specific objectives of this research are:

1) To formulate representative nonlinear boundary value problems arising in real-world applications.

2) To discretize these problems using appropriate finite difference schemes.

3) To solve the resulting nonlinear algebraic systems numerically using the computational software, Maple.

4) To compare the numerical solutions with known analytical or benchmark solutions.

5) To assess the limitations of the method and suggest improvements or alternatives.

6) To demonstrate the practical utility of the method in modeling nonlinear physical systems.

Through this work, we aim to contribute to the advancement of reliable and computationally efficient techniques for solving nonlinear BVPs and to support their broader application in scientific computing and engineering practice.

This paper is organized as follows: Section 2 reviews the literature; Section 3 outlines assumptions and limitations; Section 4 describes the study; Section 5 presents the finite difference scheme; Section 6 discusses accuracy analysis; Section 7 provides numerical examples; Section 8 examines the rate of convergence; Section 9 discusses the results; and Section 10 concludes the paper.

2. Literature Review

The numerical analysis of nonlinear boundary value problems (BVPs) arising from ordinary differential equations (ODEs) has become a crucial area of research in mathematics and applied sciences. These problems frequently emerge in real-world systems such as fluid dynamics, heat transfer, chemical reactions, population models, elastic structures, and other engineering applications, where the governing models are nonlinear and defined over specific boundary conditions. However, solving these BVPs analytically is often challenging or even impossible due to their inherent nonlinearity and complexity. To overcome such challenges, researchers have developed and employed a variety of numerical methods aimed at obtaining approximate yet accurate solutions. One of the most widely used techniques is the finite difference method (FDM), which converts differential equations into systems of algebraic equations using discrete approximations. Due to its conceptual simplicity, computational efficiency, and strong adaptability, the FDM has become a popular choice for solving linear and nonlinear BVPs.

Over the past several decades, a significant body of literature has been dedicated to the development and improvement of numerical techniques for BVPs. In particular, the fourth-order and higher-order nonlinear BVPs, which frequently arise in modeling beam deflection, viscous flows, and biological processes, have received considerable attention. For instance, Noor and Mohyud-Din applied the variational iteration method (VIM) to effectively handle fourth-order boundary value problems [1]. Hashim employed the Adomian decomposition method (ADM) to tackle both linear and nonlinear fourth-order integro-differential equations [2]. Similarly, Pandey used the finite difference method to develop numerical solutions for sixth-order linear BVPs [3]. Salama and Mansour proposed a finite difference scheme for third-order nonlinear BVPs [4], while Wang and Guo designed a compact finite difference method with a non-isotropic mesh to solve two-dimensional fourth-order nonlinear elliptic BVPs [5]. In another approach, Dang and Luan introduced an iterative technique to solve nonlinear fourth-order boundary problems [6]. Momani and Moadi utilized both the classical and modified forms of the Adomian decomposition method to obtain numerical solutions of BVPs with two-point boundary conditions [7]. Mohyud-Din and Noor investigated the homotopy perturbation method (HPM) for similar problems [8], while Liang and Jeffrey developed the homotopy analysis method (HAM) to find series solutions of fourth-order two-point boundary problems [9]. Xu contributed to the field by solving fourth-order BVPs using the variational iteration method [10], and Hayani and Casasus addressed similar problems through an approximate analytical approach using decomposition methods [11]. Additional researchers have applied spline methods, collocation techniques, and shooting methods to solve nonlinear boundary value problems effectively [12]-[15]. These studies demonstrate the continuing evolution of numerical tools and the increasing focus on improving the accuracy, stability, and computational efficiency of solutions.

In summary, the existing literature clearly highlights the wide range of techniques developed for solving BVPs and emphasizes the importance of efficient numerical approaches for tackling nonlinear ODEs. The finite difference method, in particular, remains a powerful and reliable tool in this domain. Therefore, this study focuses on using an efficient version of the finite difference method to solve nonlinear BVPs modeled by ordinary differential equations, aiming to obtain accurate numerical results while ensuring computational efficiency and robustness.

3. Assumptions and Limitations

3.1. Assumptions

1) The nonlinear boundary value problems under consideration are assumed to be well-posed, meaning they have solutions that are unique and continuously depend on the boundary value.

2) The differential equations and boundary conditions are sufficiently smooth to allow for discretization and numerical approximation using finite differences.

3) The resulting system of algebraic equations is solved using by a Maple software.

4) It is assumed that rounding and truncation errors introduced during numerical computation are within acceptable limits.

3.2. Limitations

1) The study is restricted to one-dimensional or simplified forms of nonlinear BVPs to make the implementation and analysis manageable.

2) Only stationary (time-independent) problems are considered; time-dependent PDEs are outside the scope of this research.

3) The numerical method may suffer from convergence issues for problems with strong nonlinearities or discontinuities.

4) The quality of the results depends on the step size used in discretization; very fine meshes increase accuracy but also computational cost.

5) Only Dirichlet or simple boundary conditions are considered; more complex boundary types such as Neumann or Robin conditions may require additional analysis.

4. Description of the Study

This study employs a quantitative research design, focusing on the application of numerical methods to solve nonlinear boundary value problems (BVPs) for ordinary differential equations (ODEs). The research aims to investigate the effectiveness and efficiency of finite difference methods (FDM), among other numerical approaches, for obtaining approximate solutions to these problems.

The study follows a structured procedure:

1) Problem Formulation: Initially, a set of representative nonlinear boundary value problems (both linear and nonlinear) is selected, modeled as ODEs with appropriate boundary conditions.

2) Discretization of the Domain: The continuous domain of the differential equations is discretized using a grid-based approach. The chosen method, finite difference, approximates the derivatives using finite differences at discrete grid points.

3) Numerical Implementation: The nonlinear boundary value problems are discretized using the finite difference method, and the resulting system of equations is solved using Maple software to obtain numerical approximations.

4) Error Analysis: The absolute error and relative error are calculated by comparing the numerical solutions to exact solutions (if available) or results from other established numerical methods.

5) Programming Tools: The numerical simulations and computations are implemented using MATLAB and Maple. These programming environments are well-suited for solving differential equations and performing numerical analysis, with built-in libraries for matrix operations and algebraic solvers.

6) Validation Techniques: To ensure the reliability of the results, numerical solutions are compared with known exact solutions (for cases where they exist) or validated using other numerical methods from the literature.

5. Finite Difference Scheme for Fourth Order Boundary Value Problem

The primary objective of this research is to investigate nonlinear ordinary differential equations subject to boundary conditions. Finite difference methods serve as powerful numerical tools for solving boundary value problems by replacing differential equations with corresponding difference equations through mathematical discretization. To demonstrate the fundamental concepts and implementation of the finite difference technique, we consider a general fourth-order non-linear boundary value problem (BVP) of the form:

$\frac{{\text{d}}^{4}y}{\text{d}{x}^4}=f\left(x,y\left(x\right),y^{\prime }\left(x\right),y^{″}\left(x\right),y^{‴}\left(x\right)\right);a\le x\le b$ (5.1)

With boundary conditions:

$y\left(a\right)={\alpha }_0$ ; $y^{\prime }\left(a\right)={\beta }_0$ ; $y\left(b\right)={\alpha }_1$ ; $y^{\prime }\left(b\right)={\beta }_1$ (5.2)

or

$y\left(a\right)={\alpha }_0$ ; $y^{″}\left(a\right)={\beta }_0$ ; $y\left(b\right)={\alpha }_2$ ; $y^{″}\left(b\right)={\beta }_2$ (5.3)

where $f$ is a continuous function on [a, b] and the parameters ${\alpha }_i$ and ${\beta }_i$ , where, $i=0,1,2$ are real constants. let us consider, $a = x_0 \le x_1\le x_3\le \cdots \le x_{n-1}\le x_n= b$ represents a regular partition of the interval [a, b], that is, $x_i=a+ih;i=0,1,2,\cdots ,n$ and $h=\frac{b-a}{n}$ . The points $x_1=a+h,x_2=x_1+h,\cdots ,x_n=x_{n-1}+h$ are called interior mesh points of the interval $a\le x\le b$ .

Also, for finding numerical approximate solutions we consider, $y_i=y\left(x_i\right)$ and, second order central difference approximations to the first four derivatives of y(x) are as follows:

$y^{\prime }\left(x\right)\approx \frac{y_{i+1}-y_{i-1}}{2h}$

$y^{″}\approx \frac{y_{i+1}-2y_i+y_{i-1}}{h^2}$

$y^{‴}\approx \frac{y_{i+2}-2y_{i+1}+2y_{i-1}-y_{i-2}}{2h^3}$

$y\text{'}\text{'}\text{'}\text{'}\approx \frac{y_{i+2}-4y_{i+1}+6y_i-4y_{i-1}+y_{i-2}}{h^4}$

In this method, the derivatives in the differential equation and boundary conditions are replaced with their corresponding finite difference approximations. The resulting system of algebraic equations is then solved using a Maple software. The solutions obtained at the discrete grid points represent the approximate values of the required solution to Equation (5.1) at the pivotal nodes.

6. Accuracy Analysis Scheme

The accuracy of the numerical solution obtained using the proposed finite difference method largely depends on the choice of the step size $h$ . As $h$ becomes smaller, the approximation generally becomes more accurate, thereby reducing the numerical error. A numerical solution $y\left(x_n\right)$ is said to converge to the exact solution $y_n$ if, for any given $\delta >0$ , there exists a sufficiently large positive integer $N$ such that $\left|y\left(x_n\right)-y_n\right|<\delta$ for all $n\ge N$ .

When an exact analytical solution is available, the numerical error is measured using the absolute error formula:

$e_r=\left|y\left(x_n\right)-y_n\right|.$

This approach provides a quantitative assessment of the method’s accuracy and demonstrates the effect of mesh refinement on the convergence of the numerical solution.

7. Numerical Examples

In this section, two nonlinear boundary value problems (BVPs) are considered to evaluate the accuracy of the proposed finite difference method. Numerical results are computed and presented in Tables 1-6, while the corresponding errors are illustrated graphically in Figures 1-6.

Example 1: We consider a fourth-order simplified forms of nonlinear beam deflection boundary value problem is given by

$y^{\left(iv\right)}+y^2{y}^{‴}+y^{″}=-{\sin }^{2}x\cdot \cos x$

subject to the boundary conditions:

$y\left(0\right)=0,y\left(1\right)=\sin \left(1\right),y^{″}\left(0\right)=0,y^{″}\left(1\right)=-\sin \left(1\right)$

The exact solution of the problem is $y\left(x\right)=\sin x$ . The approximate numerical results and the associated errors are computed using the proposed finite difference method and are presented in Tables 1-3. Graphical representations of the numerical solution and error behavior are shown in Figures 1-3.

Table 1. Numerical results for step size $h=0.1$ .

Values of x	Exact Solution	Numerical Solution	Absolute Error	Relative Error (%)
0.1	0.09983341665	0.09795969559	0.187372106e−2	1.876847576
0.2	0.1986693308	0.1950699157	0.35994151e−2	1.811761828
0.3	0.2955202067	0.2904877281	0.50324786e−2	1.702921995
0.4	0.3894183423	0.3833794515	0.60388908e−2	1.550746368
0.5	0.4794255386	0.4729198838	0.65056548e−2	1.356968763
0.6	0.5646424734	0.5582886701	0.63538033e−2	1.125279022
0.7	0.6442176872	0.6386646358	0.55530514e−2	0.8619836913
0.8	0.7173560909	0.7132190523	0.41370386e−2	0.5767064158
0.9	0.7833269096	0.7811088751	0.22180345e−2	0.2831556624

Table 2. Numerical results for step size $h=0.05$ .

Values of x	Exact Solution	Numerical Solution	Absolute Error	Relative Error (%)
0.1	0.09983341665	0.09795584209	0.187757456e−2	1.880707506
0.2	0.1986693308	0.1950632166	0.36061142e−2	1.815133813
0.3	0.2955202067	0.2904797130	0.50404937e−2	1.705634195
0.4	0.3894183423	0.3833717268	0.60466155e−2	1.552730019
0.5	0.4794255386	0.4729137642	0.65117744e−2	1.358245207
0.6	0.5646424734	0.5582849144	0.63575590e−2	1.125944168
0.7	0.6442176872	0 .6386633057	0.55543815e−2	0.8621901588
0.8	0.7173560909	0.7132195041	0.41365868e−2	0.5766434345
0.9	0.7833269096	0.7811098828	0.22170268e−2	0.2830270188

Table 3. Numerical results for step size $h=0.025$ .

Values of x	Exact Solution	Numerical Solution	Absolute Error	Relative Error (%)
0.1	0.09983341665	0.09983425936	8.42710000E−07	8.44116157E−04
0.2	0.1986693308	0.1986709614	1.63060000E−06	8.20760806E−04
0.3	0.2955202067	0.2955225166	2.30990000E−06	7.81638598E−04
0.4	0.3894183423	0.3894211707	2.82840000E−06	7.26314016E−04
0.5	0.4794255386	0.4794286752	3.13660000E−06	6.54241326E−04
0.6	0.5646424734	0.5646456609	3.18750000E−06	5.64516513E−04
0.7	0.6442176872	0.6442206248	2.93760000E−06	4.55994931E−04
0.8	0.7173560909	0.7173584375	2.34660000E−06	3.27117875E−04
0.9	0.7833269096	0.7833282879	1.37830000E−06	1.75954634E−04

Figure 1. Comparison of exact and numerical solutions.

Figure 2. Comparison of absolute errors for different step sizes.

Figure 3. Comparison of relative errors (%) for different step size.

Example 2: We consider a fourth-order nonlinear boundary value problem for an exponential growth in biological model is given by

$y^{\left(iv\right)}+y\cdot y^{″}+y^2=16{\text{e}}^{2x}+5{\text{e}}^{4x}$

subject to the boundary conditions:

$y\left(0\right)=1,y\left(1\right)={\text{e}}^2,y^{″}\left(0\right)=4,y^{″}\left(1\right)=4{\text{e}}^2$

The exact solution of this problem is $y\left(x\right)={\text{e}}^{2x}$ . The approximate numerical results and corresponding errors are computed using the proposed finite difference method. The outcomes are presented in Tables 4-6, and the graphical illustrations of the solution and errors are shown in Figures 4-6.

Table 4. Numerical results for step size $h=0.1$ .

Values of x	Exact Solution	Numerical Solution	Absolute Error	Relative Error (%)
0.1	1.221402758	1.221761296	0.358538e−3	0.2935460868e−1
0.2	1.491824698	1.492468218	0.643520e−3	0.4313643559e−1
0.3	1.822118800	1.823019279	0.900479e−3	0.4941933534e−1
0.4	2.225540928	2.226696784	0.1155856e−2	0.5193595793e−1
0.5	2.718281828	2.719693034	0.1411206e−2	0.5191536747e−1
0.6	3.320116923	3.321753038	0.1411206e−2	0.4250470790e−1
0.7	4.055199967	4.056959712	0.1759745e−2	0.4339477743e−1
0.8	4.953032424	4.954693538	0.1661114e−2	0.3353731326e−1
0.9	6.049647464	6.050806089	0.1158625e−2	0.1915194244e−1

Table 5. Numerical results for step size $h=0.05$ .

Values of x	Exact Solution	Numerical Solution	Absolute Error	Relative Error (%)
0.1	1.221402758	1.221494266	0.91508e−4	0.7492041376e−2
0.2	1.491824698	1.491989170	0.164472e−3	0.1102488786e−1
0.3	1.822118800	1.822349033	0.230233e−3	0.1263545494e−1
0.4	2.225540928	2.225836254	0.295326e−3	0.1326985257e−1
0.5	2.718281828	2.718641859	0.360031e−3	0.1324480031e−1
0.6	3.320116923	3.320533517	0.416594e−3	0.1254757015e−1
0.7	4.055199967	4.055647093	0.447126e−3	0.1102599141e−1
0.8	4.953032424	4.953453606	0.421182e−3	0.8503517925e−2
0.9	6.049647464	6.049940651	0.293187e−3	0.4846348514e−2

Table 6. Numerical results for step size $h=0.025$ .

Values of x	Exact Solution	Numerical Solution	Absolute Error	Relative Error (%)
0.1	1.221402758	1.2205299837	8.72774300E−04	7.14567160E−02
0.2	1.491824698	1.4911096461	7.15051900E−04	4.79313622E−02
0.3	1.822118800	1.8213228028	7.95997200E−04	4.36852526E−02
0.4	2.225540928	2.2248593792	6.81548800E−04	3.06239616E−02
0.5	2.718281828	2.7175128204	7.69007600E−04	2.82902086E−02
0.6	3.320116923	3.3193584258	7.58497200E−04	2.28454966E−02
0.7	4.055199967	4.0546528928	5.47074200E−04	1.34906837E−02
0.8	4.953032424	4.9526132169	4.19207100E−04	8.46364538E−03
0.9	6.049647464	6.0488979234	7.49540600E−04	1.23898228E−02

Figure 4. Comparison of exact and numerical solutions.

Figure 5. Comparison of absolute errors for different step sizes.

Figure 6. Comparison of relative errors (%) for different step sizes.

8. Rate of Convergence

To evaluate the efficiency and accuracy of the proposed numerical scheme, we performed a convergence analysis based on two test problems using three different step sizes: $h=0.1$ , $h=0.05$ , and $h=0.025$ . The aim is to measure how the error decreases as the step size becomes smaller. This is quantitatively expressed through the order of convergence, calculated using the standard logarithmic formula:

$p=\frac{\log \left(\frac{e_{k+1}}{e_k}\right)}{\log \left(\frac{h_{k+1}}{h_k}\right)}$

where $e_k$ and $e_{k+1}$ are the maximum errors (absolute or relative) corresponding to step sizes $h_k$ and $h_{k+1}$ , respectively. This formula provides the empirical convergence rate, indicating how closely the method approaches its theoretical accuracy.

For Problem 1, the method showed no observable convergence between $h=0.1\to 0.05$ both absolute and relative errors stagnated, with convergence orders near zero. However, a significant drop in error was observed between $h=0.05\to 0.025$ , resulting in extremely high convergence orders $p\approx 10.996$ for absolute error and $p\approx 11.12$ for relative error indicating rapid convergence once the asymptotic regime is reached. In contrast, Problem 2 exhibited smoother and more predictable behavior. Absolute errors decreased steadily, yielding a moderate convergence order of $p\approx 1.98$ . However, when the step size was further refined from $h=0.05\to 0.025$ the convergence order dropped to $p\approx -0.9649$ indicating a breakdown in convergence likely due to numerical instability or stiffness in the problem. In terms of relative error, the method initially exhibited near-quadratic convergence with $p\approx 1.97$ between $h=0.1\to 0.05$ , but the order significantly declined to $p\approx -2.4289$ for $h=0.05\to 0.025$ , reflecting numerical saturation or divergence at finer grids. Despite these challenges, the method retains acceptable accuracy at moderate resolutions, and its effectiveness may be enhanced with stabilization techniques when applied to stiff or highly nonlinear problems.

Table 7. Table of convergence behavior.

Problem	Error Type	Step Interval	Order p	Behavior
1	Absolute	$h=0.1\to 0.05$	–0.0014	No convergence
1	Absolute	$h=0.05\to 0.025$	10.996	High convergence
1	Relative	$h=0.1\to 0.05$	–0.0030	No convergence
1	Relative	$h=0.05\to 0.025$	11.12	Rapid convergence
2	Absolute	$h=0.1\to 0.05$	1.98	Moderate convergence
2	Absolute	$h=0.05\to 0.025$	–0.9649	No convergence
2	Relative	$h=0.1\to 0.05$	1.97	Near quadratic convergence
2	Relative	$h=0.05\to 0.025$	–2.4289	No significant convergence at the finer grid

9. Discussion of Results

The results of the numerical solution of the nonlinear boundary value problem (BVP) using the proposed finite difference method are presented in Tables 1-3 and Tables 4-6, with corresponding graphical illustrations shown in Figures 1-3, and Figures 4-6. The convergence behavior associated with these results is summarized in Table 7. Numerical solutions were computed for step sizes $h=0.10$ , $h=0.05$ and $h=0.025$ , and were compared against the exact solution to assess the performance and accuracy of the method. The graphical representations clearly demonstrate that the numerical solutions closely approximate the exact solution, with improved accuracy observed for the smaller step size $h=0.025$ . As the step size decreases, the numerical values align more closely with the exact values, indicating a higher level of precision and better convergence. This is a typical characteristic of well-behaved numerical schemes. To quantitatively evaluate the accuracy of the method, both absolute error and relative error have been calculated. The absolute error, defined as the direct difference between the exact and numerical solutions, shows a noticeable reduction as the step size becomes smaller. This suggests that finer discretization leads to better local approximations at each node.

In addition, the relative error, which provides a percentage-based measure of accuracy, has been analyzed and illustrated in Figure 3 and Figure 6. The relative error values consistently decrease across the solution domain when the smaller step size $h=0.025$ is used. However, while the proposed finite difference method generally demonstrates a high rate of convergence and accuracy as the step size decreases, certain cases exhibit stagnation or divergence, particularly at finer grids. This observation suggests that although the method performs well for a wide range of nonlinear boundary value problems, its accuracy and stability may be affected by factors such as problem stiffness, round-off errors, or numerical saturation when the grid becomes excessively refined. Overall, the results confirm that the proposed method is effective, with convergence towards the exact solution validated through both tabular data and graphical evidence. The method’s performance supports its practical applicability to nonlinear BVPs, while also highlighting the importance of addressing potential limitations in stiff or highly nonlinear cases.

10. Conclusions

In this study, we have explored the numerical solution of fourth-order nonlinear boundary value problems (BVPs) using the finite difference method (FDM). The primary objective was to develop an efficient and reliable numerical scheme capable of approximating solutions to nonlinear BVPs, particularly in situations where analytical methods are intractable due to the complexity of the equations.

The proposed finite difference method was successfully formulated and implemented to discretize the governing differential equations. Numerical solutions were obtained for various step sizes, and their accuracy was assessed through comparison with known exact solutions. Both absolute and relative errors were evaluated to quantify deviations, while the analysis of convergence behavior demonstrated that the accuracy of the solution improves as the step size decreases. The graphical and tabular comparisons further confirmed the stability, consistency, and reliability of the approach. Overall, the results validate the effectiveness of the finite difference method as a practical and efficient computational tool for solving nonlinear BVPs. Its simplicity, adaptability, and low computational overhead make it highly suitable for a wide range of applications in mathematical modeling, engineering, and the physical sciences.

As a direction for future research, this method can be extended to more complex and higher-dimensional nonlinear problems, including time-dependent systems. Investigating adaptive mesh strategies, hybridizing FDM with other numerical schemes (e.g., finite element or spectral methods), and applying the method to stiff or singularly perturbed BVPs are promising areas. Moreover, the integration of parallel computing can significantly enhance performance for large-scale simulations. These advancements will further strengthen the role of FDM as a robust tool in computational mathematics.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Noor, M.A. and Mohyud-Din, S.T. (2008) Variational Iteration Method for Solving Fourth-Order Boundary Value Problems. Applied Mathematics and Computation, 195, 591-596.
[2]	Hashim, I. (2006) Adomian Decomposition Method for Solving Fourth-Order Integro-Differential Equations. Applied Mathematics and Computation, 181, 324-330.
[3]	Pandey, R.K. (2001) Numerical Solution of Sixth-Order Boundary Value Problems Using Finite Difference Method. Journal of Computational and Applied Mathematics, 134, 303-316.
[4]	Salama, M.M.A. and Mansour, S.M. (2008) Finite Difference Methods for Solving Third-Order Nonlinear Boundary-Value Problems. Applied Mathematics and Computation, 204, 395-404.
[5]	Wang, B. and Guo, Z. (2006) A Compact Finite Difference Method with Non-Isotropic Mesh for Solving Fourth-Order Nonlinear Elliptic Problems. Computers & Mathematics with Applications, 51, 949-962.
[6]	Dang, H.S. and Luan, N.T. (2007) An Iterative Method for Solving Fourth-Order Nonlinear Boundary Value Problems. Applied Mathematics and Computation, 189, 1951-1961.
[7]	Momani, S. and Moadi, M. (2007) Adomian Decomposition Method for Solving Fourth-Order Boundary Value Problems. Journal of Computational and Applied Mathematics, 207, 185-193.
[8]	Mohyud-Din, S.T. and Noor, M.A. (2009) Homotopy Perturbation Method for Solving Fourth-Order Boundary Value Problems. Communications in Nonlinear Science and Numerical Simulation, 14, 1901-1908.
[9]	Liang, X. and Jeffrey, D.J. (2006) Homotopy Analysis Method for Solving Fourth-Order Boundary Value Problems. Applied Mathematics and Computation, 172, 500-510.
[10]	Xu, L. (2007) Variational Iteration Method for Solving Fourth-Order Boundary-Value Problems. International Journal of Computer Mathematics, 84, 1347-1356.
[11]	Hayani, S. and Casasus, R. (2005) Decomposition Method for Solving Fourth-Order Boundary Value Problems. Applied Mathematics Letters, 18, 917-922.
[12]	Bhatti, M.I. and Bracken, P. (2007) Spline Solutions of Nonlinear Boundary Value Problems in Mathematical Modeling. Mathematical and Computer Modelling, 45, 148-160.
[13]	Dehghan, M. (2005) Numerical Solution of a Class of Nonlinear Boundary Value Problems Using Collocation Methods. Mathematical and Computer Modelling, 41, 1395-1407.
[14]	Ascher, U.M., Mattheij, R.M.M. and Russell, R.D. (1995) Numerical Solution of Boundary Value Problems for Ordinary Differential Equations. SIAM. https://doi.org/10.1137/1.9781611971231
[15]	Keller, H.B. (1968) Numerical Methods for Two-Point Boundary-Value Problems. Dover Publications.