Comparison of the method of variation of parameters to semi-analytical methods for solving nonlinear boundary value problems in engineering

Abstract Solutions to nonlinear boundary value problems modelling physical phenomena in engineering applications have traditionally been approximated using numerical methods. More recently, several semi-analytical methods have been developed and used extensively in diverse engineering applications. This work compares the method of variation of parameters to semi-analytical methods for solving nonlinear boundary value problems arising in engineering. The accuracy and efficiency of the method of variation of parameters are compared to those of two widely used semi-analytical methods, the Adomian decomposition method and the differential transformation method, for three practical engineering boundary value problems: (1) the deflection of a cantilevered beam with a concentrated load, (2) an adiabatic tubular chemical reactor which processes an irreversible exothermal reaction, and (3) the electrohydrodynamic flow of a fluid in an ion drag configuration in a circular cylinder conduit. The accuracy and convergence of each method is investigated using the error remainder function. The method of variation of parameters is significantly more efficient than both the semi-analytical methods and traditional numerical methods while maintaining comparable accuracy. Unlike the semi-analytical methods, the efficiency of variation of parameters is independent of the nonlinearity. Variation of parameters is shown to be an attractive alternative to semi-analytical methods and traditional numerical methods for solving boundary value problems encountered in engineering applications in which solution efficiency is important.


Introduction
Mathematical modelling of physical phenomena in all branches of science and engineering frequently results in boundary value problems governed by nonlinear differential equations. These problems generally do not have closed-form, exact analytical solutions. However, because of their many important applications, significant efforts have been and continue to be made to accurately approximate solutions to these problems. Early methods developed to approximate solutions to nonlinear boundary value problems (BVPs) were strictly numerical. Whereas analytical methods are symbolic in nature, numerical techniques involve reformulating problems such that their solutions can be approximated with arithmetic operations whose error can be estimated [1]. Commonly used numerical methods for solving BVPs include shooting methods, Runge-Kutta methods, finite-difference methods, and integral equation methods. Extensive treatment of these and other numerical methods for solving BVPs are widely available [2,3].
More recently, several so-called "semi-analytical" methods for solving differential equations have been developed. Among these are the variational iteration method (VIM), the homotopy perturbation method (HPM), the homotopy analysis method (HAM), the Adomian decomposition method (ADM), and the differential transformation method (DTM). These methods have been used extensively in recent years and several texts and countless articles have been devoted to investigating their application to solve differential equations in numerous engineering applications from diverse fields such as solid mechanics, heat transfer, fluid mechanics, dynamics, and biomedicine [3][4][5][6][7][8][9][10]. The advantages of many of these semianalytical approaches over numerical methods are their direct applicability to both linear and nonlinear equations without requiring linearization, discretization, or perturbation [11]. Additionally, they can be used to prove existence of solutions. However, these methods result in series solutions which must be truncated and can cause convergence problems. Although typically more accurate, semi-analytical methods are generally slower than numerical methods. This can be problematic in cases where solution efficiency is critical. An important example is inverse problems, in which the governing equation must be solved numerous times for optimization purposes [12].
A potential efficient and accurate alternative to semianalytical and traditional numerical methods for solving nonlinear BVPs is the method of variation of parameters (VPM). This method was developed by Lagrange to solve linear, nonhomogeneous ordinary differential equations [13]. However, it can also be implemented to solve nonlinear differential equations in which the nonlinearity is found in the inhomogeneity [14][15][16][17]. The primary contribution of the present work is the comparison of the accuracy and efficiency of VPM to common semi-analytical methods. For reference, these methods will also be compared to traditional finite difference methods. It will be shown that VPM is a versatile and practical alternative for efficiently and accurately solving nonlinear BVPs that arise from models of physical phenomena in a wide variety of engineering applications.
Following an overview of how VPM, ADM, and DTM are implemented to solve nonlinear BVPs, solutions to four BVPs are presented using these methods in addition to traditional finite-difference methods. The first is a simple example with a known analytical solution so that the accuracy of each method can be quantified. Then, the solutions to three recently investigated engineering applications of nonlinear problems are shown. These include modelling (1) the deflection of a cantilevered beam with a concentrated load, (2) an adiabatic tubular chemical reactor which processes an irreversible exothermal reaction, and (3) the electrohydrodynamic flow of a fluid in an ion drag configuration in a circular cylinder conduit. These problems represent diverse fields of engineering whose solutions using semi-analytical methods have recently been published. Their governing equations have distinct linear and nonlinear portions which is used to demonstrate the versatility of the method of variation of parameters. The accuracy and efficiency of each method are compared for the three applications. Accuracy and convergence of each method is investigated using the error remainder function. The method of variation of parameters is shown to be an attractive alternative to semi-analytical methods and traditional numerical methods for solving boundary value problems encountered in engineering applications, especially for problems in which solution efficiency is an important consideration.

Method of variation of parameters
This section will demonstrate how the method of variation of parameters can be used to solve nonlinear differential equations in which the nonlinearity is in the inhomogeneity. For the sake of convenience and simplicity, the following description of the method is limited to second order differential equations. However, this method can be extended to higher order equations [18]. Consider the following equation where f is a nonlinear function of the independent variable x as well as the dependent variable y and/or its derivative y ′ . The solution to this equation consists of the sum of the complementary and particular solutions. The complementary solution is the solution to the homogeneous equation corresponding to (1) that may be found using traditional ordinary, linear differential equation solution techniques [19]. By the principle of superposition, the complementary solution is yc = c 1 y 1 + c 2 y 2 where c 1 and c 2 are arbitrary constants and y 1 and y 2 form a fundamental set of solutions of the homogeneous equation. Herein lies the restriction that the homogeneous equation must be linear. Therefore, this method can be used for any BVP of the form . . y (n−1) )︁ [18]. The VPM method assumes that there exist functions u 1 and u 2 such that a particular solution has the form yp = u 1 y 1 + u 2 y 2 [19]. The unknown functions u 1 and u 2 are found by substituting the assumed particular solution into (1) while imposing the condition u ′ 1 y 1 + u ′ 2 y 2 = 0, which results in a system of equations with unknowns u ′ 1 and u ′ 2 . This system of equations can be solved to find expressions for u ′ 1 and u ′ 2 which can, in turn, be separated and integrated to yield the desired functions.
The solution to (1) is therefore given by where the constants c 1 and c 2 are determined from the boundary conditions. Equation (4) is an exact solution to (1). However, due to the implicit nature of this solution, an iterative approach is required to approximate y. Numerical integration must be used during the iteration process and, if f is a function of y ′ , finite difference methods are required to approximate the derivative. The accuracy of this method is therefore limited by the accuracy of the numerical methods employed. Details regarding the comparisons between exact solutions and those obtained by VPM are available in [20].

Adomian decomposition method
The Adomian decomposition method (ADM) is a common semi-analytical method used to solve nonlinear ordinary and partial differential equations [21]. It is applied by dividing the equation into linear and nonlinear parts and applying the inverse of the highest order linear derivative operator to both parts. The boundary conditions along with the source term are identified as the zeroth component of the series solution while the remaining terms are found by decomposing the nonlinear part into a series in which the terms are calculated recursively using Adomian polynomials. For the sake of consistency, the following description of this method is limited to second order equations of the form (1). However, the approach is applicable to nth order equations. Suppose that (1) is expressed as where L = d 2 dx 2 is the second-order derivative operator, R = p (x) d dx + q (x) is the linear operator for the remainder of the linear portion of the equation, and N is the nonlinear operator N (y) = f (︁ x, y, y ′ )︁ . It is assumed that the inverse operator L −1 exists and can be taken as a definite integral of a function as follows.
Moving R (y) to the right side of (5) and applying the inverse operator to both sides yields where y 0 = y (x0) + y x 0 is the first component of the decomposition series given by The nonlinear function N (y) is represented by the infinite series Nn (9) where the components Nn are the Adomian polynomials which are defined as [21] Nn (y0, y 1 , . . . , yn) The remaining components of the series in (8) are determined recursively from the equation ADM is a versatile method in that it can be used for linear and nonlinear ordinary and partial differential equations. However, it requires significant computational time to calculate the Adomian polynomials for higher order approximations. Additionally, this method is rapidly convergent in only very small regions while the convergence rate is much slower over wider regions [22].

Differential transformation method
The differential transformation method is another frequently used semi-analytical method. This method was first proposed by Zhou [23] to solve linear and nonlinear initial value problems arising in electric circuit analysis and has subsequently been used to solve numerous equations in a wide variety of applications [7,24]. DTM is based on the Taylor series expansion and, like ADM, results in a series solution to the differential equation. Consider again the differential equation (1). An approximate solution to this equation is given by the series where the coefficients Y (0), Y (1). . . Y (N) are obtained by applying the differential transform (13) to both sides of the differential equation. The result is an algebraic equation whose unknowns are the coefficients of the series solution (12). Common operation properties of the one-dimensional differential transform can be found in [23,25]. Equation (1) can be rearranged as follows and the differential transform operations applied, resulting in the following recursive relation where Y (k), P (k), Q (k), and F (k) are the differential transformations of y (x), p (x), q (x), and the nonlinear function f (︁ x, y, y ′ )︁ , respectively. This equation can be used along with the transformed boundary conditions to find the values of Y (k) required in the series solution. As is the case for ADM, DTM results in a truncated series solution which only provides an accurate approximation of the true solution in a very small region.

Example BVP: Comparison to analytical solution
The following example will demonstrate the use of all three methods to solve BVPs. Consider the following boundary value problem.
The analytical solution is If (16) is rearranged as it can be treated as an equation of the form of (1) with p (x) = q (x) = 0 and f = e x − y and can be solved using numerical and semi-analytical methods as presented below.

Solution by the method of variation of parameters
The complementary solution to (18) is yc = c 1 x + c 2 where y 1 = x and y 2 = 1. Using the method of variation of parameters as outlined in Section 2 and applying the boundary conditions to find the constants c 1 and c 2 results in the following solution to (18).
This solution is implicit and must be solved iteratively and the integrals must be evaluated numerically.

Solution by the differential transformation method
Application of the differential transform operators to (18), where the differential transform of e x is 1/k!, gives the following recurrence relation.
From the differential transform given by (13), the transformed boundary conditions are Utilizing the recurrence relation in (23) and the first transformed boundary condition, Y (k) for k ≥ 2 are found and the following series solution to (18) is obtained from (12).
The second transformed boundary condition is used to solve for the unknown Y (1) = y ′ (0).

Results and discussion
In order to properly compare the accuracy and efficiency of the method of variation of parameters to the semianalytical methods, all three methods were implemented in MATLAB [26], a high-level scripting language. Equation (19) was solved iteratively with an initial guess of y (x) = 0 for 101 discrete values of x equally spaced over the domain. Simpson's rules were used to perform the numerical integration. Convergence was considered achieved when the Euclidean norm of the difference between the current and previous iteration values of y fell below 1 × 10 −9 . Further decreasing this criterion had no effect on the solution for all 15 significant figures of the double precision values. Convergence was achieved after 25 iterations which took 0.0338 seconds. The average percent error between the VPM solution and the exact solution at the discrete values was 3.89 × 10 −7 %. To solve (18) in MATLAB using the Adomian decomposition method, the terms in the series solutions were solved for symbolically in terms of the unknown y ′ 0 . At each iteration, the boundary condition y ′ (1) = 0 was enforced to calculate a value for y ′ 0 using a numeric solver and the series approximation of y was calculated over the domain at 101 discrete values of x. Convergence was considered achieved when the Euclidean norm of the difference between the current and previous iteration values of yfell below 1 × 10 −12 . Further decreasing this criterion had no effect on the solution. Convergence was achieved in 3.24 seconds with 9 terms in the series solution. The average percent error between the ADM solution and the exact solution at the discrete values was 3.54 × 10 −14 %. Implementation of the differential transformation method to solve (18) in MATLAB followed the same procedure as that used for ADM. The terms in the series were solved for in terms of the unknown constant Y (1), which was approximated at each iteration by applying the second boundary condition in (24). At each iteration, the series approximation of y was calculated over the domain at 101 discrete values of x. Convergence was considered achieved when the Euclidean norm of the difference between the current and previous iteration values of y fell below 1 × 10 −16 . Further decreasing this criterion had no effect on the solution. Convergence was achieved in 0.047 seconds with 21 terms in the series solution. The average percent error between the DTM solution and the exact solution at the discrete values was 3.14 × 10 −14 %.
In order to compare these methods to traditional numerical methods, (18) was also solved with the "bvp4c" routine provided in MATLAB. This routine is a finite difference method (FDM) that implements the three-stage Lobatto IIIa collocation formula to solve boundary value problems. The solution to (18) was approximated at 101 discrete values of x. The relative tolerance was set to 1 × 10 −7 below which there was no change to the solution. The computational time required was 0.278 seconds and the average percent error between the FDM solution and the exact solution at the discrete values was 1.35 × 10 −7 %. Table 1 shows a comparison of the computational time required to achieve convergence and the average percent error for each of the methods. It also compares the results of each method at six discrete points with 10 significant figures. Figure 1 shows a plot of the solutions at discrete points. All methods produced highly accurate results relative to the exact solution. The accuracies of the two semianalytical methods were roughly the same and were much greater than those of the method of variation of parameters and the finite difference method, which also had comparable accuracies. VPM was the most efficient of all methods. The required computational time was much lower than ADM and FDM and slightly lower than DTM. The accuracy of VPM and FDM can be increased by using smaller discretization, which would decrease the efficiency. The accuracy of VPM can also be increased by using more accurate numerical quadrature methods. The following practical engineering applications show that the method of variation of parameters is, in general, significantly faster than semi-analytical methods and traditional numerical methods for soling boundary value problems while also maintaining a high level of accuracy.

Application 1: Deflection of a cantilevered beam with a concentrated load
Bisshopp and Drucker [27] obtained the following nonlinear boundary value problem as a model of the deflection of a cantilever beam under a concentrated load, where θ is the angle (radians) of the slope of the deflection curve at a point along the beam, s is the nondimensional arc length of the beam, and λ is the ratio of the load to the flexural rigidity of the beam (m −2 ). This model is applicable to beams that experience large deflections and is more realistic than those typically used in texts on mechanics of materials [28] because it does not neglect the square of the first derivative in the curvature formula when applying the Bernoulli-Euler theorem. Bisshopp and Drucker solved the equation using elliptical integrals [27]. Na obtained numerical solutions using the finite difference method, the method of reduced physical parameters, and the method of invariant embedding [2]. More recently, a numerical solution was obtained by the collocation method with Haar wavelets [29]. The exact solution to (26) is unknown.

Solution by the method of variation of parameters
The complementary solution to (26) is θc = c 1 s + c 2 where θ 1 = s and θ 2 = 1. Using the method of variation of parameters as outlined in Section 2 and applying the boundary conditions to find the constants c 1 and c 2 results in the following solution to (26), where f (s, θ) = −λs cos θ is the nonlinear function.

Solution by the differential transformation method
Taking the differential transform of both sides of (26) results in the following recursive relation where G (k) is the differential transform of cos θ which can be found in [30]. The boundary conditions are transformed as follows The approximate series solution is obtained by substituting the expressions calculated from (30) into (12) and applying the first transformed boundary condition, resulting in an equation in terms of the unknown θ (0). This unknown is approximated from the second transformed boundary condition.

Results and discussion
Equation (26) was solved in MATLAB using each of the methods with λ = 8, consistent with [2,27,29]. The same procedures as those described in Section 5.4 were used to implement the four solution methods. Convergence for VPM was considered achieved when the Euclidean norm of the difference between the current and previous iteration values of θ fell below 1 × 10 −13 . Convergence was achieved after 103 iterations which took 0.0889 seconds. The convergence criterion for ADM was 1 × 10 −5 . The Adomian polynomials were calculated with the aid of the MAT-LAB function "allVl1" [31]. Convergence was achieved in 24.7 seconds with 11 terms in the series solution. Due to the extremely slow rate of convergence, the convergence criterion for DTM was set to 1×10 −4 . Convergence was achieved in 20.8 minutes with 28 terms in the series solution. Increasing the convergence criterion decreased the required computational time but also decreased the accuracy of the solution. The relative tolerance of the "bvp4c" routine for FDM was 1 × 10 −9 . The computational time required was 0.312 seconds. Table 2 shows a comparison of the computational time required to achieve convergence for each of the methods and compares the results at six discrete points with 10 significant figures. Figure 2 shows a plot of the solutions at discrete points.
Because the exact solution to (26) is unknown, the error remainder function is used as a measure of how well the approximate solutions satisfy the governing equation. The error remainder function for (26) is where θn is the approximate solution to (26) after n iterations. The average error remainder function over the interval 0 ≤ s ≤ 1 at different values of n is a measure of the accuracy of the approximate solution as well as the rate of convergence. Figure 3 compares the average error remainder function as a function of iteration for the VPM, ADM, and DTM methods. Finite difference methods were employed to approximate θ ′′ n in (32). The accuracy and rate of convergence for VPM is comparable to that of DTM while ADM has better accuracy and convergence.

Application 2: Adiabatic tubular chemical reactor
A mathematical model of the temperature distribution in an irreversible exothermic chemical reaction processed in an adiabatic tubular chemical reactor under steady state conditions is given by [32] where y is the unknown nondimensional steady state temperature, x is the dimensionless length, λ is the Péclet number which represents the ratio of the rate of advection to the rate of diffusion, µ is the Damköhler number which represents the ratio of the chemical reaction rate to the rate of convective mass transport, and β is the dimensionless adiabatic temperature rise. The existence of the solution to this problem over various ranges has been investigated [33] and approximate solutions have been obtained using Adomian's method for Hammerstein integral equations and the shooting method [34], the Sinc-Galerkin method [35], and the Sinc-collocation method [36].

Solution by the method of variation of parameters
The linear, homogeneous equation corresponding to (33) is solved by finding the roots of the auxiliary equation m 2 − λm = 0 , resulting in the complementary solution yc = c 1 + c 2 exp (λx) where y 1 = 1 and y 2 = exp (λx).
Using the method of variation of parameters as outlined in Section 2 and applying the boundary conditions to find the constants c 1 and c 2 results in the following solution to (33), where f (y) = λµ (y − β) exp (y) is the nonlinear inhomogeneity.

Solution by the Adomian decomposition method
Application of the inverse operator to both sides of (33) and using (7) yields where L −1 is given by (6) with x 0 = 0. The first boundary condition combined with the first two terms on the right side of (35) provide the first term in the series solution y 0 = y (0) + λx y (0). The recursive relation is therefore where Nn (y) = Nn (y0, y 1 , . . . , yn) are the Adomian polynomials of the nonlinear operator N (y) = (y − β) exp (y). The series solution can be found in terms of the unknown y (0), which is found by applying the second boundary condition to the approximate solution.

Solution by the differential transformation method
Taking the differential transform of both sides of (33) results in the following recursive relation where W (k) is the differential transform of exp (y) which can be found in [30]. The boundary conditions are transformed as follows The approximate series solution is obtained by substituting the expressions calculated from (37) into (12) and applying the first transformed boundary condition, resulting in an equation in terms of the unknown y (0). This unknown is approximated from the second transformed boundary condition.

Results and discussion
Equation (33) was solved in MATLAB using each of the methods with λ = 5, µ = 0.05, and β = 0.53, consistent with [34][35][36]. The same procedures as those described in Section 5.4 were used to implement the four solution methods. The convergence criterion for VPM was 1 × 10 −12 which was achieved in 0.0284 seconds after 8 iterations.
Due to the slow rate of convergence of ADM, the convergence criterion set to 1 × 10 −5 . Convergence was achieved in 25.7 minutes with 6 terms in the series solution. The Adomian polynomials were calculated with the aid of the MAT-LAB function "allVl1" [31]. Increasing the convergence criterion decreased the required computational time but also decreased the accuracy of the solution. The convergence criterion for DTM was 1 × 10 −7 which was achieved in 25.2 seconds with 20 terms in the series solution. An underrelaxation factor of 0.01 was required to achieve convergence. The relative tolerance of the "bvp4c" routine for FDM was 1 × 10 −6 . The computational time required was 0.225 seconds. Table 3 shows a comparison of the computational time required to achieve convergence for each of the methods and compares the results at six discrete points with 10 significant figures. Figure 4 shows a plot of the solutions at discrete points. The error remainder function for (33) is where yn is the approximate solution to (33) after n iterations. Figure 5 compares the average error remainder function over the interval 0 ≤ x ≤ 1 as a function of iteration for the VPM, ADM, and DTM methods. The accuracy and rate of convergence for ADM and DTM are comparable. The initial accuracy of VPM is better than the semi-analytical methods but converges after few iterations. This may be a result of the limitation of using finite-difference methods to approximate y ′′ n and y ′ n in (39) for VPM.

Application 3: Electrohydrodynamic flow
Electrohydrodynamics is the study of the dynamics of electrically charged fluids. Mathematical modelling of electrohydrodynamic flow is important in the design and analysis of flowmeters, accelerators, pumps, and magnetohydrodynamic generators. The electrohydrodynamic flow of a fluid in an "ion drag" configuration in a circular cylindrical conduit is governed by the nonlinear differential equation [37] where w is the nondimensional fully-developed fluid velocity, r is the nondimensional radial distance from the center of the cylindrical conduit (the radius of the tube has been scaled to one), H is the Hartmann number which is the ratio of electromagnetic force to viscous force, and α is a dimensionless parameter related to the pressure gradient, the ion mobility, and the current density at the inlet of the conduit. Approximate solutions to (40) have been obtained using traditional numerical methods (Runge-Kutta/shooting method) [37] as well as with semianalytical methods such as DTM [38] and the homotopy perturbation method [39]. The existence and uniqueness of the solution to (40) is examined in [40].

Solution by the method of variation of parameters
The linear, homogeneous equation corresponding to (40) is an Euler-Cauchy equation which can be solved by finding the roots of the auxiliary equation m 2 = 0 , resulting in the complementary solution wc = c 1 + c 2 ln (r) where w 1 = 1 and w 2 = ln (r). Using the method of variation of parameters as outlined in Section 2 and applying the boundary conditions to find the constants c 1 and c 2 results in the following solution to (40), )︀ is the nonlinear inhomogeneity.

Solution by the Adomian decomposition method
In (40) it is more efficient to keep the (︀ 1/r )︀ dw/dr term on the left side of the equation such that the equation can be written as where the operator L is defined by )︂ (43) and the inverse operator is therefore Applying the inverse operator to both sides of (42) and enforcing the boundary conditions leads to the following recursive relation where Nn (w) = Nn (w0, w 1 , . . . , wn) are the Adomian polynomials of the nonlinear operator N (w) = −H 2 (︀ 1 − w/ (1 − αw) )︀ . By using the operator in (43) instead of the second-order derivative operator, application of the boundary conditions results in w 0 = 0 thereby eliminating the need to use a solver to find w 0 . Equation (45) is used with (8) to determine the approximate series solution.

Solution by the differential transformation method
Taking the differential transform of both sides of (40) results in the following recursive relation where F (k) is the differential transform of 1 − w/ (1 − αw) which is given by The boundary conditions are transformed as follows The approximate series solution to (40) is obtained by substituting the expressions calculated from (46) into (12) and applying the first transformed boundary condition, resulting in an equation in terms of the unknown w (0). This unknown is approximated from the second transformed boundary condition.

Results and discussion
Equation (40) was solved in MATLAB using each of the methods with H = 0.25 and α = 0.25, consistent with [37][38][39]. The same procedures as those described in Section 5.4 were used to implement the four solution methods. The convergence criterion for VPM was 1 × 10 −13 which was achieved in 0.0281 seconds after 8 iterations. The convergence criterion for ADM was 1 × 10 −6 which was achieved in 1.57 seconds with 4 terms in the series solution. The Adomian polynomials were calculated with the aid of the MAT-LAB function "allVl1" [31]. The convergence criterion for DTM was 1 × 10 −13 which was achieved in 0.398 seconds with 9 terms in the series solution. The relative tolerance of the "bvp4c" routine for FDM was 1 × 10 −8 . The computational time required was 0.381 seconds. Table 4 shows a comparison of the computational time required to achieve convergence for each of the methods and compares the results at six discrete points with 10 significant figures. Figure 6 shows a plot of the solutions at discrete points.
where wn is the approximate solution to (40) after n iterations. Figure 7 compares the average error remainder function over the interval 0 ≤ r ≤ 1 as a function of iteration for the VPM, ADM, and DTM methods. The accuracy and rate  of convergence for ADM and DTM are comparable. The initial accuracy of VPM comparable to that of the other methods but converges quickly while the accuracy of ADM and DTM continue to increase as the number of iterations increases. Again, this may be a result of the limitation of using finite-difference methods to approximate w ′′ n and w ′ n in (49) for VPM.

Conclusion
This work has demonstrated the method of variation of parameters to be an accurate and very efficient method for solving nonlinear boundary value problems arising in engineering applications. The accuracy and efficiency of this method were compared to those of two semi-analytical methods, the Adomian decomposition method and the dif-ferential transformation method, as well as to traditional finite difference methods. Each of these methods were used to approximate solutions to BVPs encountered in diverse engineering applications whose governing equations had distinct linear and nonlinear portions. Generally, the accuracy of the semi-analytical methods was slightly better than that of the method of variation of parameters. At worst, the solutions from the different methods were equivalent up to at least five significant figures. However, the computational efficiency of the method of variation of parameters was greater than that of all other methods and was generally independent of the governing equation. In contrast, the efficiency of the semi-analytical methods was highly dependent on the nonlinearity of the BVP. The accuracy of the method of variation of parameters can be increased by using smaller discretization of the domain and by employing more accurate numerical quadrature methods to approximate the integrals. The method of variation of parameters is an attractive alternative to semi-analytical methods and traditional numerical methods for solving boundary value problems encountered in engineering applications, especially for problems in which solution efficiency is an important consideration.