J. Nig. Soc. Phys. Sci. 5 (2023) 1087 Journal of the Nigerian Society of Physical Sciences Numerical Solution of Second Order Fuzzy Ordinary Differential Equations using Two-Step Block Method with Third and Fourth Derivatives Kashif Hussain, Oluwaseun Adeyeye∗, Nazihah Ahmad School of Quantitative Sciences, Universiti Utara Malaysia, Kedah, Malaysia Abstract Fuzzy differential equation models are suitable where uncertainty exists for real-world phenomena. Numerical techniques are used to provide an approximate solution to these models in the absence of an exact solution. However, existing studies that have developed numerical techniques for solving second-order fuzzy ordinary differential equations (FODEs) possess an absolute error accuracy that could be improved. Therefore, this article developed a more accurate higher derivative self-starting block scheme for the numerical solution of second-order FODEs with fuzzy initial and boundary conditions imposed. Linear block approach using Taylor series expansion is adopted for the derivation of the proposed method and the basic properties are established using the definitions of stability and consistency for block methods. According to the numerical results, when compared to the exact solution in terms of absolute error, the new method proposed in this article outperformed existing numerical methods. It is thus concluded that the proposed method is effective for solving second-order FODEs directly. DOI:10.46481/jnsps.2023.1087 Keywords: Fuzzy initial value problem, Fuzzy boundary value problem, Second order, Two-step, Block method, Linear, Nonlinear Article History : Received: 24 September 2022 Received in revised form: 20 January 2023 Accepted for publication: 12 February 2023 Published: 04 April 2023 © 2023 The Author(s). Published by the Nigerian Society of Physical Sciences under the terms of the Creative Commons Attribution 4.0 International license (https://creativecommons.org/licenses/by/4.0). Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. Communicated by: B. J. Falaye 1. Introduction Second-order differential equations have many applications, especially in the field of engineering, biology, chemistry, elec- tronics, physics, etc. Unfortunately, unpredictable scenarios may be encountered which introduced the concept of uncer- tainty [1] and the application of fuzzy derivatives in fuzzy dif- ferential equations (FDEs) to handle these situations [2]. There are three differentiations used to describe the differential or deriva- tive of a fuzzy function. The first is the Hukuhara derivative ∗Corresponding author tel. no: +60 49286354 Email address: adeyeye@uum.edu.my (Oluwaseun Adeyeye) (H-derivative), which was introduced in [3], the second is the Seikkala derivative introduced in [4], and the third is the gen- eralized derivative (g-derivative) introduced in [5]. This study focuses on the H-derivative in order to define the differential equations considered in this article, which follows the defini- tion by the authors whose results were considered for compari- son in the numerical examples with the newly developed block method. The second-order FODE of the form given in the equation below is considered in this article, ŷ′′(x) = f (x, ŷ(x), ŷ′(x)),∀x ∈ [a, b] (1) 1 Hussain et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1087 2 From Equation 1, ŷ′′(x) = d 2 ŷ dt2 = f (x, ŷ(x), ŷ ′(x)) is a H- derivative and ŷ is a fuzzy function of crisp variable x. Since the function is fuzzy, there exist solutions known as lower and upper solutions because the parametric form of the α-level is given as ŷ ′′ (x,α) = f (x, ŷ(x,α), ŷ′(x,α)),∀α ∈ [0, 1], where f = min { f (x, ŷ(x,α), ŷ′(x,α)) } and f = max { f (x, ŷ(x,α), ŷ ′ (x,α)) } . The above types of problems in the parametric form of fuzzy function may be difficult to solve directly, and sometimes it is not possible to obtain exact solutions. As a result, researchers were interested in employing various numerical approaches to obtain an approximate solution for second-order FODEs. Sev- eral types of numerical methods developed by numerous re- searchers for second-order FODEs with initial and boundary conditions include the homotopy analysis method in [6, 7], de- composition method [8], Laplace and differential transforma- tion method in [9, 10], least-square method [11], and Runge- Kutta method in [12-14]. The biggest drawback of these ap- proaches is the reduction of the second-order FODEs to the sys- tem of first-order FODEs, which leads to computational burden and also impacts solution accuracy. To bypass the rigor of re- duction, block methods were introduced for the direct solution of second-order FODEs in [15-17]. However, due to the order of the block methods developed by these studies, it is observed that there is still room to improve the accuracy of their obtained results in terms of absolute error. Hence, the motivation of this study is to develop a new block method with the presence of two higher derivative terms with the aim of obtaining better accuracy. In comparison to existing methods, the newly de- veloped method has the advantages of better accuracy, being self-starting, and incurring a low computational burden in the development and implementation of the block method. The following is how this article is structured: The essential definitions for fuzzy set theory are presented in Section 2, and the construction of the two-step block method with third and fourth derivatives is presented in Section 3 with the use of the linear block approach. Section 4 highlights the block method’s properties, Section 5 considers linear and nonlinear numerical examples, and Section 6 concludes the article. 2. Preliminaries This section recalls some definitions which will be adopted in this article. The section discusses basic definitions of trian- gular fuzzy numbers, trapezoidal fuzzy numbers, fuzzy set sup- port, α-level set, and Hukuhara differential. These concepts are required to establish the different parameters of the crisp the- ory’s uncertain behavior. These concepts play an important role when fuzzy differential equations model real-life situations. Definition 1: Triangular Fuzzy Number [18] Consider three numbers (µ, v, w) ∈ R3,µ ≤ v ≤ w, then M(x) denotes the triangular fuzzy number given as: M(x,µ, v, w) =  0, x < µ x−µ v−µ , µ ≤ x ≤ v w−x w−v , v < x ≤ w 0, x > w (2) The corresponding α-level set is defined as Mα = [ µ + α (v −µ) , w −α(w − v) ] ,α ∈ [0, 1]. (3) Definition 2: Trapezoidal Fuzzy Numbers [18] Consider four numbers (µ, v, w,δ) ∈ R4,µ ≤ v ≤ w ≤ δ, then the trapezoidal fuzzy number M(x) is given as: M(x,µ, v, w,δ) =  0, x < µ x−µ v−µ , µ ≤ x < v 1, v ≤ x ≤ w w−x w−v , w < x ≤ δ 0, x > δ (4) The corresponding α-level set is defined as Mα = [ µ + α (v −µ) ,δ−α(δ− w) ] ,α ∈ [0, 1]. (5) Definition 3: Fuzzy Set Support [18] A set  has fuzzy set support with X universal set defined as, S upp(Â) = { x ∈ X|MÂ(x) > 0 } (6) It contains all elements in X which have membership degree of fuzzy element greater than zero. Definition 4: α-Level Set [18] Consider that, M ∈ R f , the α-level set is defined as, Mα = {x ∈ R|M(x) > 0} , α ∈ [0, 1]cl(suppM), α = 0 , (7) with its closed, bounded interval [M(x), M(x)]. M(x) and M(x) are lower and upper bound of Mα respectively. Definition 5: Hukuhara Differential [3] A function f : (u, v) → R f is called H-differentiable, if for h > 0 sufficiently small, then H-difference f (x) − f (x − h), f (x + h) − f (x) exists and ∃ an element f ′(x) ∈ R f such that, lim h→0 f (x) − f (x − h) h = lim h→0 f (x + h) − f (x) h = f ′(x). (8) Then f ′(x) is called the H-derivative of f at x. 2 Hussain et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1087 3 3. Methodology Given that the second-order FODE defined in Equation 1 be a mapping f : R f → R f and ŷ0 ∈ R f with α-level set ŷ0 ∈ (̂ y(0,α), ŷ(0,α) )α α ,α ∈ [0, 1]. The partition of the has the set of grid points 0 = x0 < x1 < x2 <,...,< xn = X with exact solution as ( Ŷ (xn,α) )α α = ( Ŷ (xn,α), Ŷ (xn,α) )α α and approxima- tion solution also denoted as (̂ y(xn,α) )α α = (̂ y(xn,α), ŷ(xn,α) )α α at which points, h = X−x0n , xn = x0 + nh, 0 ≤ n ≤ N. The two-step linear block method with the presence of third and fourth derivatives in second-order form is stated below as, (̂ yn+η )α α =  1∑ v=0 (ηh)v v! ŷ(v)n + 2∑ d=0  2∑ v=0 ψdvη f (d) n+v   α α ,η = 1, 2 (9) with the first derivative expression for the block method form given as (̂ y′n+η )α α = ̂y′n + 2∑ d=0  2∑ v=0 ωdvη f (d) n+v   α α ,η = 1, 2 (10) Expanding Equations 9 and 10 produces the expressions in Equations 11, 12, 13, and , 14. (̂ yn+1 )α α =  ŷn + ĥy ′ n + [ψ001 fn + ψ011 fn+1 + ψ021 fn+2 +ψ101 f ′ n + ψ111 f ′ n+1 + ψ121 f ′ n+2 + ψ201 f ′′ n +ψ211 f ′′ n+1 + ψ221 f ′′ n+2]  α α (11) (̂ yn+2 )α α =  ŷn + 2ĥy ′ n + [ψ002 fn + ψ012 fn+1 + ψ022 fn+2 +ψ102 f ′ n+1 + ψ112 f ′ n+1 + ψ122 f ′ n+2 + ψ202 f ′′ n +ψ212 f ′′ n+1 + ψ222 f ′′ n+2]  α α (12) (̂ y ′ n+1 )α α =  ŷ′n + [ω001 fn + ω011 fn+1 + ω021 fn+2 + ω101 f ′ n +ω111 f ′ n+1 + ω121 f ′ n+2 + ω201 f ′′ n + ω211 f ′′ n+1 +ω221 f ′′ n+2]  α α (13) (̂ y ′ n+2 )α α =  ŷ′n + [ω002 fn + ω012 fn+1 + ω022 fn+2 + ω102 f ′ n +ω112 f ′ n+1 + ω122 f ′ n+2 + ω202 f ′′ n + ω212 f ′′ n+1 +ω222 f ′′ n+2]  α α (14) By applying Taylor series expansions (̂ y(x + h; α) )α α =  n∑ i=0 hi i! f i(x; α) α α (15) which is given in [19] to expand each term in Equations 11-14 yields (̂ yn+ j )α α = (̂ y(xn + jh; α) )α α =  n∑ i=0 ( jh)i i! f i(xn; α) α α , j = 0, 1, 2, (16) (̂ yn+ j )α α =  ŷ(xn; α) + jĥy ′(xn; α) + ( jh)2 2! ŷ′′(xn; α) + ( jh)3 3! ŷ′′′(xn; α) + .... + ( jh)n n! ŷn(xn; α)  α α . (17) After that, the unknown coefficients ψdvn and ωdvn are ob- tained from ψdvn = A−1 B and ωdvn = A−1 D, where A =  1 1 1 0 0 0 0 0 0 0 h 2h 1 1 1 0 0 0 0 h 2 2! 22 h2 2! 0 h 2h 1 1 1 0 h 3 3! 23 h3 3! 0 h2 2! 22 h2 2! 0 h 2h 0 h 4 4! 24 h4 4! 0 h3 3! 23 h3 31 0 h2 2! 22 h2 2! 0 h 5 5! 25 h5 5! 0 h4 4! 24 h4 4! 0 h3 3! 23 h3 3! 0 h 6 6! 26 h6 6! 0 h5 5! 25 h5 5! 0 h4 4! 24 h4 4! 0 h 7 7! 27 h7 7! 0 h6 6! 26 h6 6! 0 h5 5! 25 h5 5! 0 h 8 8! 28 h8 8! 0 h7 7! 27 h7 7! 0 h6 61 26 h6 6!  α α , B = α (ηh)2 2! (ηh)3 3! (ηh)4 4! (ηh)5 5! (ηh)6 6! (ηh)7 7! (ηh)8 8! (ηh)9 9! (ηh)10 10!  α , D =  ηh (ηh)2 2! (ηh)3 3! (ηh)4 4! (ηh)5 5! (ηh)6 6! (ηh)7 7! (ηh)8 8! (ηh)9 9!  α α . Therefore,  ψ001 ψ011 ψ021 ψ101 ψ111 ψ121 ψ201 ψ211 ψ221  α α =  19h2 60 h2 5 −h2 60 911h3 20160 −16h3 315 113h3 20160 53h4 20160 h4 80 −11h4 20160  ,  ψ002 ψ012 ψ022 ψ102 ψ112 ψ122 ψ202 ψ212 ψ222  α α =  76h2 105 128h2 105 2h2 35 34h3 315 −32h3 315 −2h3 315 2h4 315 16h4 315 0  ,  ω001 ω011 ω021 ω101 ω111 ω121 ω201 ω211 ω221  α α =  5669h 13440 64h 105 −42h 13440 303h2 4480 −1h2 8 47h2 4480 169h3 40320 8h3 315 −41h3 40320  ,  ω002 ω012 ω022 ω102 ω112 ω122 ω202 ω212 ω222  α α =  41h 105 128h 105 41h 105 2h2 35 0 −2h2 35 1h3 315 16h3 315 1h3 315  . The obtained values of the coefficients are substituted in Equations 11-14 which is the required two-step block method 3 Hussain et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1087 4 with the presence of third and fourth derivatives as given below.(̂ yn+1 )α α = ŷn + ĥy ′ n + h 2 [ 19 60 fn + 1 5 fn+1 − 1 60 fn+2 ] +h3 [ 911 20160 gn − 16 315 gn+1 + 113 20160 gn+2 ] +h4 [ 53 20160 mn + 1 80 mn+1 − 11 20160 mn+2 ] , (̂ yn+2 )α α = ŷn + 2ĥy ′ n + h 2 [ 76 105 fn + 128 105 fn+1 + 2 35 fn+2 ] +h3 [ 34 315 gn − 32 315 gn+1 − 2 315 gn+2 ] + h4 [ 2 315 mn + 16 315 mn+1 ] , (18) (̂ y′n+1 )α α = ŷ ′ n + h [ 5669 13440 fn + 64 105 fn+1 − 421 13440 fn+2 ] +h2 [ 303 4480 gn − 1 8 gn+1 + 47 4480 gn+2 ] +h3 [ 169 40320 mn + 8 315 mn+1 − 41 40320 mn+2 ] , (̂ y′n+2 )α α = ŷ ′ n + h [ 41 105 fn + 128 105 fn+1 + 41 105 fn+2 ] +h2 [ 2 35 gn − 2 35 gn+2 ] + h3 [ 1 315 mn + 16 315 mn+1 + 1 315 mn+2 ] (19) where g = d f (x,α)d x , m = d2 f (x,α) d x . The block method in Equation 18 has corrector form,( A0Ŷn+k )α α = ( A1Ŷn−k )α α + h ( B1Ŷ ′ n−k )α α + h2 ( C0 Fn+k + C 1 Fn−k )α α +h3 ( D0Gn+k + D 1Gn−k )α α + h4 ( E0 Mn+k + E 1 Mn−k )α α where, A0 = ( 1 0 0 1 )α α , A1 = ( 0 1 0 1 )α α , B1 = ( 0 1 0 2 )α α , C0 = ( 1 5 −1 60 128 105 2 35 )α α , C1 = ( 0 1960 0 76105 )α α , D0 = ( −16 315 113 20160 −32 315 −2 315 )α α , D1 = ( 0 91120160 0 34315 )α α , E0 = ( 1 80 −11 20160 16 315 −2 315 )α α , E1 = ( 0 5320160 0 2315 )α α , Ŷn+k = (̂ yn+1 ŷn+2 )α α , Ŷn−k = (̂ yn−1 ŷn )α α , Ŷ′n−k = (̂ y′n−1 ŷ′n )α α , Fn+k = ( fn+1 fn+2 )α α , Fn−k = ( fn−1 fn )α α , Gn+k = ( gn+1 gn+2 )α α Gn−k = ( gn−1 gn )α α , Mn+k = ( mn+1 mn+2 )α α , Mn−k = ( mn−1 mn )α α . 4. Properties of the Proposed Method This section will first mention the required definitions and theorems to investigate the properties of the developed two-step third-fourth derivative scheme, and thereafter apply these theo- rems and definitions to the method. 4.1. Convergence and Stability Properties Theorem 1: A block method is convergent iff it is consistent and zero- stable. [22] Proof The aim of the proof is to show that zero stability and con- sistency are necessary conditions for convergence. Suppose that the block method defined in Equation 9 is convergent, the first condition for zero-stability follows by considering Equation 1 with a trivial solution ŷ(x) = 0. Applying Equation 9 to this problem yields the difference equation̂yn+η − 1∑ v=0 (ηh)v v! ŷ(v)n − 2∑ d=0  2∑ v=0 ψdvη f (d) n+v   α α ,η = 1, 2 (20) Since the method is assumed to be convergent, for any x > 0, then lim h→0 nh→0 ŷn+η = 0 (21) for all solutions of Equation 20 satisfying ŷs = ςs(h), s = 0, 1, ..., k − 1 where lim h→0 ŷs = 0 (22) Let ψ = reiφ be a root of the first characteristic polynomial P(ψ) = 0, r ≥ 0, 0 ≤ φ ≤ 2π. It can be verified then that the numbers ŷn+η = hr n cos(nφ) (23) define a solution to Equation 20 satisfying Equation 22. If φ = 0, φ , π, then ŷn+η − ŷn − ŷ ′ n sin2φ = h2r2n (24) Since the left-hand side of this identity converges to 0 as h → 0, n →∞, nh = x the same must be true of the right-hand side; therefore, lim n→∞ ( x n )∞ r2n = 0 (25) This implies that r ≤ 1. In other words, it is proven that any root of the first characteristic polynomial of (9) lies in the closed unit disc. Note that any root of the first characteristic polynomial of Equation 9 that lies on the unit circle must be simple. For the other condition, which is consistency, let us first show that C0 = 0. Consider Equation 1 with trivial solution, ŷ(x) = 1. Applying Equation 9 to this problem yields the differ- ence equation Equation 20. Choose ŷs = 1, s = 0, 1, ..., k − 1. Given that by hypothesis the method is convergent, it is deduced that lim h→0 ŷs = 1 (26) Since in the present case ŷn is independent of the choice of h, Equation 26 is equivalent to saying that lim h→∞ ŷn = 1, (27) 4 Hussain et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1087 5 and passing to the limit n → ∞ in Equation 20, it is deduced that αk + αk−1+, ..., +α0 = 0. (28) Recalling the definition of C0, Equation 28 is equivalent to C0 = 0 (i.e. P(1) = 0). To show that C1 = 0, consider Equation 1 with trivial solu- tion, ŷ(x) = x. Applying Equation 9 to this problem yields the difference equation in Equation 20. For a convergent method every solution of Equation 20 satisfying lim h→0 ςs(h) = 0, s = 0, 1, ..., k − 1 (29) where ŷs = ςs(h), s = 0, 1, ..., k − 1, must also satisfy lim h→0 ŷn+η = x. (30) Since according to the previous theorem zero-stability is necessary for convergence, we may take it for granted that the first characteristic polynomial P(ψ) of the method does not have multiple roots on the unit circle |ψ| = 1, therefore P ′ (1) = kαk+, ..., +2α2 + α1 , 0. (31) Let the sequence (xn) N n = 0 be defined by ŷn = Knh, where K = ψdkη + ψd(k−1)η+, ..., +ψd2η + ψd1η + ψd0η kαk+, ..., +2α2 + α1 . (32) This sequence clearly satisfies Equation 30 and is the solu- tion of Equation 20. Furthermore, Equation 31 implies that x = ŷ(x) = lim h→0 nh=x ŷn+η = lim h→0 nh=x Knh = K x (33) C1 = (kαk+, ..., +2α2 + α1) −(ψdkη + ψd(k−1)η+, ..., +ψd2η + ψd1η + ψd0η) = 0. (34) Equivalently, P ′ (1) = σ(1). Thus, since the necessary con- ditions in terms of zero-stability and consistency is satisfied, so the block method is convergent. Definition 6: Consistency [20] A block method is consistent if it has order ρ ≥ 1. Definition 7: Zero-Stability [20] A block method with matrix difference equation in the fol- lowing form A0Ŷn+k = A 1Ŷn−k + B 1Ŷ ′′ n−k + B 2Ŷ ′′ n−k + · · · + B 1Ŷ (m−1)n−k +hm ( C0Ŷ mn+k + C 1Ŷ mn−k )α α + h(m+1) ( D0Ŷ (m+1)n+k + D 1Ŷ (m+1)n−k )α α +h(m+2) ( E0Ŷ (m+2)n+k + E 1Ŷ (m+2)n−k )α α , (35) with Ŷ an+k = (̂ yan+1, ŷ a n+2, ..., ŷ a n+k )T and Ŷ an−k = (̂ yan1(k−1), ŷ a n−(k−2), ..., ŷ a n )T , is zero-stable if the first char- acteristic polynomial takes form P(ψ) = det(ψv A 0 − A1), (36) and the root of P(ψ) = 0 satisfy |ψv| ≤ 1, v = 1, ..., , k. Definition 8: Region of Absolute Stability [26] To obtain the polynomial for the absolute stability region of the block method. The expressions for the corrector take the form: det  −(w)k + A1 + q  k∑ j=0 B jwk− j  + q2  k∑ j=0 theC jwk− j  +q3  k∑ j=0 D jwk− j  + q4  k∑ j=0 E jwk− j    α α , q = λh. The absolute stability region is then obtained by plotting the polynomial roots using the boundary locus technique. If the obtained roots of the polynomial lie in the unit circle, then the block method is absolutely stable and its region is called the region of absolute stability. Note that large absolute stability regions mean that large time-step size can be used during the implementation of the method to solve the differential equation [27-29]. Definition 9: A-stable According to [20], a numerical method is said to be A-stable if its region of absolute stability contains the whole of the left- hand half-plane. Definition 10: L-stable According to [20] a general linear multistep method is L- stable if it is A-stable and, in addition, when applied to the scalar test equation ŷ ′ = λy, λ is a complex constant with Reλ < 0, it yields ŷn+1 = R(hλ)̂yn, where, |R(hλ)|→ 0 as Re(hλ) →∞. However, a clause is encountered as given in the following def- inition Definition 11 According to [21] an A-stable linear multistep method can- not have an order greater than two. Therefore, based on Definition 8, the properties of A-stability and L-stability cannot be explored for the block methods devel- oped in this article. This is because the block method devel- oped have order greater than two. Hence, the stability property with respect to choosing a stepsize value is limited to just ab- solute stability alone. Although, much attention was not placed on choosing h-values from the stability region because the h- values were chosen the same as the authors for comparison. These definitions for block methods in crisp form is adopted to the proposed method for FODEs to prove the convergence properties for the proposed method in the next subsection. 4.2. Convergence and Stability Analysis of Proposed Method Order and Error constant The linear operator associated with Equation 9 is defined as: L(̂y(x), h) = ̂yn+η − 1∑ v=0 (ηh)v v! ŷ(v)n + 2∑ d=0  2∑ v=0 ψdvη f (d) n+v   α α , η = 1, 2, (37) 5 Hussain et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1087 6 with L(̂y(x), h) = C0̂y(xn) + C1ĥy′(xn) + C2h2̂y′′(xn) + ... +Cz+1h z+1̂yz+1(xn) + Cz+2h z+2̂yz+2(xn)  α α . The method is said to be of order z if C0 = C1 = · · · = Cz = Cz+1 = 0, Cz+2 , 0, and Cz+2 is the error constant. Following the approach by [28], The order of the two-step third-fourth derivatives block method with corrector Equation 18 is nine with an error constant C11 = (3.8174e − 08, 7.617e − 08) T , and the order of the deriva- tive part ten with an error constant C12 = (6.5076e − 08, −7.6349e − 09) T . The derivative formu- lae will be used to obtain the first derivative term in Equation 1. Expressing the corrector scheme 18 as blocks using previous definitions for the block methods. A simple iteration has been implemented to approximate the value of ŷn+1 and ŷn+2. In the code, we iterate the corrector to convergent and the convergence test employed, and the order of the correctors in nine [23] Zero-Stability Applying above definition in fuzzy form for the proposed method gives P(ψ) = det(ψv A 0 − A1)αα, (38) P(ψ) = ∣∣∣∣∣∣ψv ( 1 0 0 1 ) − ( 0 1 0 1 )∣∣∣∣∣∣ α α . The root of P(ψ) = 0 satisfies the condition |ψv| ≤ 1, v = 1, 2. Convergence The proposed method is convergent because it is zero stable and consistent. Absolute Stability Region The polynomial of the proposed block method to plot its region of absolute stability is obtained as:  ∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣ ( w 0 0 w2 ) + ( 0 1 0 1 ) + q ( 0 1 0 2 ) +q2 [( 1w 5 1w2 60 128w 105 2w2 35 ) + ( 0 1960 0 76105 )] +q3 [( −16w 315 113w2 20160 −32w 315 −2w2 315 ) + ( 0 91120160 0 34315 )] +q4 [( 1w 80 −11w2 20160 16w 315 0 ) + ( 0 5320160 0 2315 )] ∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣  α α , (39) R(w) =   11q8 396900 − 37q7 88200 + 19q6 52920 + 44q5 4725 + 13q4 3600 − 2q3 35 − 9q2 35 + 1]w 3 + [ 43q 8 793800 + 53q7 52920 + 508327q6 33868800 − 31547q5 793800 − 35639q4 100800 − 44q3 45 − 20597q2 20160 − 2q − 1  w  α α The absolute stability region is thus plotted as shown in Fig- ure 1, which implies that large time-stepsizes can be utilised with the method. From Figure 1, it is seen that for the absolute stability region, all the roots of polynomial lie on the unit circle. Figure 1. Absolute stability region of proposed method 5. Results This section details the application of the developed block method for the solution of second-order (linear and nonlinear) FODEs (FIVPs and FBVPs) and the obtained results are com- pared with the exact solution and existing methods. Compar- isons between exact and approximate solutions are shown in tables and graphs. x−axis shows the value of the approximation solution, y−axis show the value of α-level values, Ŷ, Ŷ are the lower and upper bounds of the exact solution re- spectively, ŷ, ŷ are the lower and upper bounds of the approximate solution respectively, E = ∣∣∣∣Ŷ − ŷ∣∣∣∣ computes the absolute error of the lower bound ap- proximation, E = ∣∣∣∣∣Ŷ − ŷ ∣∣∣∣∣ computes the absolute error of the upper bound ap- proximation, h is the step size, TSBM: Two-step Block Method with Third and Fourth Deriva- tives, EBHDEF: Extended Block Hybrid Backward Differentiation Formula [16], BDF: Block Differentiation Formula [15], BBDF: Block Backward Differentiation Formula [15], OOMB: Optimization of One-Step Block Method [17], RK5: Runge Kutta Method Order Five [14], OHAM: Optimal Homotopy Asymptotic Method [7], FDM: Finite Difference Method [30]. Example 1. Given the second-order linear FIVP ŷ′′(x) = −̂y(x), ŷ(0,α) = 0, ŷ′(0,α) = (0.9 + 0.1α, 1.1 − 0.1α), with exact solution Ŷ (x,α) = (0.9 + 0.1α) sin(x), Ŷ (x,α) = (1.1 − 0.1α) sin(x), 6 Hussain et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1087 7 Table 1. Lower and Upper solution of Example 1 α TSBM E EBHDEF E BBDF E BDF E h = 0.1 h = 0.01 h = 0.01 h = 0.01 0 0.0000e+00 2.8094e-11 5.4048e-08 3.0991e-08 0.2 0.0000e+00 2.8719e-11 5.5249e-08 3.1647e-08 0.4 0.0000e+00 2.9343e-11 5.6450e-08 3.2335e-08 0.6 0.0000e+00 2.9966e-11 5.7651e-08 3.3023e-08 0.8 0.0000e+00 3.0592e-11 5.8853e-08 3.3711e-08 1 0.0000e+00 3.1216e-11 6.0054e-08 3.4399e-08 α TSBM E EBHDEF E BBDF E BDF E h = 0.1 h = 0.01 h = 0.01 h = 0.01 0 1.1102e-16 3.4337e-11 5.4048e-08 3.7838e-08 0.2 1.1102e-16 3.3713e-11 5.5249e-08 3.7151e-08 0.4 1.1102e-16 3.3089e-11 5.6450e-08 3.6463e-08 0.6 0.0000e+00 3.2464e-11 5.7651e-08 3.5775e-08 0.8 0.0000e+00 3.1840e-11 6.1255e-08 3.5087e-08 1 0.0000e+00 3.1216e-11 6.0054e-08 3.4399e-08 Figure 2. Numerical solution of Example 1 with Lower/Upper solution and at x = 1, Ŷ (1,α) = [ Y (1,α), Y (1,α) ] , 0 ≤ α ≤ 1. The results obtained for Example 1 are shown in Table 1 and Figure 2 displays the complete iterations graph with stepsize h = 0.1 and h = 0.01 partition of the time interval x ∈ [0, 1]. It is observed from Table 1 that the approximate solution obtained by new proposed method in comparison to the exact solution in terms of absolute error is very impressive, as it give same results as the exact solution at certain points. The results are graphically shown in Figure 2. In the figure the behaviour of the linear FIVP solution is seen to monotonically increase as shown in the graph. This follows from the property that a func- tion’s output will not appear more than once during the course of a monotonically rising interval. It is worth noting that y(x) rises in lockstep with x. The exact and approximate solutions are also compared using the graph and it shows the approxi- mate solution completely overlapping the exact solution which indicates high accuracy of the proposed method. Table 2. Lower and Upper solution of Example 2 α TSBM E EBHDEF E BBDF E BDF E h = 0.1 h = 0.01 h = 0.01 h = 0.01 0 6.661338e-16 9.8449e-14 2.4250e-10 1.5988e-10 0.2 3.663736e-15 3.4927e-13 5.7971e-10 3.9122e-10 0.4 1.532108e-14 9.7144e-13 1.2016e-09 3.7933e-09 0.6 5.129230e-14 2.2859e-12 2.2597e-09 2.6125e-09 0.8 1.498801e-13 6.4525e-12 3.9207e-09 6.6967e-08 1 3.850253e-13 4.7628e-12 6.3971e-09 1.1110e-08 α TSBM E EBHDEF E BBDF E BDF E h = 0.1 h = 0.01 h = 0.01 h = 0.01 0 1.345191e-13 4.2267e-11 4.07084e-08 1.26440e-07 0.2 7.431833e-12 2.6623e-11 2.98235e-08 9.1889e-08 0.4 3.907541e-12 1.5982e-11 2.13238e-08 7.34552e-08 0.6 2.774669e-12 9.0485e-11 1.48046e-08 3.52946e-08 0.8 9.001688e-13 8.1274e-11 9.93257e-09 1.51728e-08 1 3.854694e-13 4.7628e-12 6.39707e-09 1.11097e-08 Figure 3. Numerical solution of Example 2 with Lower/Upper solution Example 2. Given the second-order non-linear FIVP ŷ′′(x) = −(̂y′(x))2, ŷ(0,α) = (α, 2 −α), ŷ′(0,α) = (1 + α, 3 −α), with exact solution Ŷ (x,α) = ln((xα + x + 1)eα), Ŷ (x,α) = ln((3x − xα + 1)eα−2), and at x = 1, Ŷ (1,α) = [ Y (1,α), Y (1,α) ] , 0 ≤ α ≤ 1. The results obtained for Example 2 are shown in Table 2 and Figure 3 displays the complete iterations graph with stepsize h = 0.1 and h = 0.01 partition of the time interval x ∈ [0, 1]. It is observed from Table 2 that the approximate solution ob- tained by the new proposed method in comparison to the exact solution in terms of absolute error is very impressive. Just as the previous example, the results graphically shown in Figure 3 are monotonically increasing showing the behaviour of the nonlinear FIVP. Likewise, the approximate solution completely 7 Hussain et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1087 8 Table 3. Lower and Upper solution of Example 3 α Exact Solution TSBM E EBHDEF E h = 0.1 h = 0.1 0 -0.100004086851013030 1.94289e-16 4.131e-07 0.2 -0.080004095094799887 8.32667e-17 4.137e-07 0.4 -0.060004103338586523 1.59594e-16 4.141e-07 0.6 -0.040004111582373492 6.93889e-17 4.149e-07 0.8 -0.020004119826159905 2.56739e-16 4.149e-07 1 -0.000004128069946763 1.35559e-16 4.161e-07 α Exact Solution TSBM E EBHDEF E h = 0.1 h = 0.1 0 0.100003908832573600 2.77555e-17 9.094e-03 0.2 0.080003917076360453 1.52655e-16 9.094e-03 0.4 0.060003925320147089 4.85722e-17 1.267e-01 0.6 0.040003933563933947 1.66533e-16 8.459e-02 0.8 0.020003941807720582 6.24500e-17 4.291e-02 1 -0.000004128069946763 1.35559e-16 4.161e-07 overlaps the exact solution which indicates high accuracy of the proposed method. Example 3. Given the second-order linear FBVP ŷ′′(x) + ŷ(x) + x = 0, ŷ(0,α) = ŷ(1,α) = (0.1α−0.1, 0.1−0.1α), with exact solution Ŷ (x,α) = −x + (0.1α− 0.1) cos(x) + (1.13376 + 0.054630α) sin(x) Ŷ (x,α) = −x + (0.1 − 0.1α) cos(x) + (1.24303 − 0.054630α) sin(x) and at x = 1, Ŷ (1,α) = [ Y (1,α), Y (1,α) ] , 0 ≤ α ≤ 1. The results obtained for Example 3 are shown in Table 3 and Figure 4 displays the complete iterations graph with stepsize h = 0.1 partition of the time interval x ∈ [0, 1]. From Table 3 and the graph in Figure 4, impressive monotono- cally dereasing results are still observed. The absolute error accuracy is high compared with the existing EBHDEF method and the overlapping behaviour of the approximate solution with the exact solution is evident. Example 4. Given the second-order non-linear FBVP ŷ′′(x) = − [̂y′(x)]2 ŷ(x) , x ∈ [0, 1], ŷ(0,α) = (0.9+0.1α, 1.1−0.1α), ŷ(1,α) = (0.9+0.1α, 2.1−0.1α), with exact solution Ŷ (x,α) = √ 1.4 + 0.1α √ 0.1(9 + α)2 14 + α + 2x, Ŷ (x,α) = √ 1.6 − 0.1α √ −0.1(−11 + α)2 −16 + α + 2x, Figure 4. Numerical solution of Example 3 with Lower/Upper solution Table 4. Lower and Upper solution of Example 4 α TSBM E EBHDEF E FDM E h = 0.1 h = 0.008 h = 0.008 0 0.000000e+00 0 0 0.2 4.4408920e-16 2.57e-06 9.27e-07 0.4 2.4424906e-15 2e-06 8.55e-07 0.6 1.5321077e-14 1.26e-06 5.92e-07 0.8 9.6207486e-12 5.88e-07 2.94e-07 1 0.000000e+00 0 0 α TSBM E EBHDEF E FDM E h = 0.1 h = 0.1 h = 0.008 0 0.000000e+00 0 0 0.2 4.4408920e-16 2.05e-06 8.15e-07 0.4 8.8817841e-15 1.63e-06 7.65e-07 0.6 1.7763568e-15 1.03e-06 5.35e-07 0.8 1.3322676e-15 4.87e-07 2.67e-07 1 0.000000e+00 0 0 and at x = 1, Ŷ (1,α) = [ Y (1,α), Y (1,α) ] , 0 ≤ α ≤ 1. The results obtained for Example 4 are shown in Table 4 and Figure 5 displays the complete iterations graph with stepsize h = 0.1 and h = 0.008 partition of the time interval x ∈ [0, 1]. It is observed from Table 4 that the approximate solution obtained by the new proposed method in comparison to the ex- act solution in terms of absolute error is very impressive as it give same results as the interval boundaries. The results are graphically shown in Figure 5 and the behaviour of the nonlin- ear FBVP solution is seen to monotonically increase. The com- parison of the exact and approximate solutions on the graph also shows high accuracy as the plots overlap. This indicates the high accuracy of the proposed method. In addition, the time in seconds required to compute the approx- imate solution of the numerical examples is given in the table below. The program code was written with MATLAB 2015a on a laptop with 8GB RAM and Intel Core i5-3427U CPU. 8 Hussain et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1087 9 Figure 5. Numerical solution of Example 4 with Lower/Upper solution Table 5. Time Taken to Compute Approximate Solutions α Example 1 Example 2 Example 3 Example 4 time/sec time/sec time/sec time/sec 0 0.4767 2.0844 0.9204 0.4457 0.2 1.3682 1.9451 1.5024 0.4097 0.4 1.3742 2.0163 1.4550 0.4072 0.6 1.0563 1.9995 1.3675 0.3741 0.8 1.8364 2.1498 1.4603 0.3982 1 0.8937 2.0415 1.4031 0.3813 α Example 1 Example 2 Example 3 Example 4 time/sec time/sec time/sec time/sec 0 0.9037 2.2057 1.3941 0.3928 0.2 0.8487 2.0648 1.3404 0.3877 0.4 0.8703 2.7846 1.4259 0.4361 0.6 0.8650 2.1738 1.4331 0.2778 0.8 0.8650 1.9303 1.4310 0.3267 1 0.8937 2.0415 1.4031 0.3813 6. Conclusion The major objective of this research to enhance the accu- racy of the solution (in terms of absolute error) by developing a numerical technique for solving second order FODEs (FIVPs and FBVPs) directly. As a result, this article developed a two- step block method for second-order FODEs with the presence of third and fourth derivatives. The proposed method outper- forms other methods discovered in the literature as shown in the tables and graphs of the numerical results obtained. In ad- dition, the method eliminates the requirement for complicated subroutines in conventional methods that require starting values or predictors. The proposed block method has proven to be a viable strategy with increased accuracy for solving both linear and nonlinear FIVPs and FBVPs. The method developed us- ing linear block approach with low computational complexity also satisfied all convergence conditions for the block methods. Hence, the proposed method in this article is more suitable for obtaining the approximate solutions of second order FIVPs and FBVPs. Acknowledgments We thank the referees for the positive enlightening com- ments and suggestions, which have greatly helped us in making improvements to this paper. References [1] S. S. Mohammed, “On fuzzy soft set-valued maps with application”, J. Nig. Soc. Phy. Sci. 2 (2020) 26. https://doi.org/10.46481/jnsps.2020.48 [2] L. A. Zadeh & S. Chang, “On fuzzy mapping and con- trol”, IEEE Trans. Syst. Man Cybern. 2 (1972) 30. https://doi.org/10.1142/9789814261302 0012 [3] M. L. Puri & D. A. Ralescu, “Differentials of fuzzy functions”, J. Math. Anal. Appl. 91 (1983) 552. https://doi.org/10.1016/0022- 247X(83)90169-5 [4] S. Seikkala, “On the fuzzy initial value problem”, Fuzzy Sets Syst. 24 (1987) 319. https://doi.org/10.1016/0165-0114(87)90030-3 [5] B. Bede & S. G. Gal, “Generalizations of the differentia- bility of fuzzy-number-valued functions with applications to fuzzy differential equations”, Fuzzy Sets Syst. 151 (2005) 581. https://doi.org/10.1016/j.fss.2004.08.001 [6] A. F. Jameel, M. Ghoreishi, & A. I. M. Ismail, “Approximate solution of high order fuzzy initial value problems”, J. Uncertain Syst. 8 (2014) 149. [7] A. Jameel, A. Shather, N. Anakira, A. Alomari, & A. Saaban, “Com- parison for the approximate solution of the second-order fuzzy nonlinear differential equation with fuzzy initial conditions”, Math. Stat. 8 (2020) 527. https://doi.org/10.13189/ms.2020.080505 [8] Ö. Akın, T. Khaniyev, Ö. Oruç, & I. Türkşen, “An algorithm for the solu- tion of second order fuzzy initial value problems”, Expert Syst. Appl. 40 (2013) 953. https://doi.org/10.1016/j.eswa.2012.05.052 [9] E. ElJaoui, S. Melliani, & L. S. Chadli, “Solving second-order fuzzy dif- ferential equations by the fuzzy Laplace transform method”, Adv. Differ. Equ. 2015 (2015) 1. https://doi.org/10.1186/s13662-015-0414-x [10] N. Salamat, M. Mustahsan, & M. M. Saad Missen, “Switch- ing point solution of second-order fuzzy differential equations using differential transformation method”, Math. 7 (2019) 231. https://doi.org/10.3390/math7030231 [11] N. J. B. Pinto, E. Esmi, V. F. Wasques, & L. C. Barros, “Least square method with quasi linearly interactive fuzzy data: fitting an HIV dataset”, Proceedings of the International Fuzzy Systems Association World Congress, Springer (2019) 177. https://doi.org/10.1007/978-3-030- 21920-8 17 [12] A. G. Fatullayev, E. Can, & C. Köroğlu, “Numerical solution of a bound- ary value problem for a second order fuzzy differential equation”, TWMS J. Pure Appl. Math. 4 (2013) 169. [13] T. Jayakumar, K. Kanagarajan, & S. Indrakumar, “Numerical solution of Nth-order fuzzy differential equation by Runge-Kutta method of order five”, Int. J. Math. Anal. 6 (2012) 2885. [14] A. Jameel, N. Anakira, A. Alomari, I. Hashim, & M. Shakhatreh, “Nu- merical solution of n-th order fuzzy initial value problems by six stages Range Kutta method of order five”, Int. J. Electr. Comput. Eng. 10 (2019) 6497. http://dx.doi.org/10.22436/jnsa.009.02.26 [15] T. K. Fook, & Z. B. Ibrahim, “Block backward differentiation formulas for solving second order fuzzy differential equations”, MATEMATIKA: Malaysian J. Ind. Appl. Math. 33 (2017) 215. https://doi.org/10.11113/matematika.v33.n2.868 [16] K. J. Audu, A. Ma’Ali, U. Mohammed, & A. Yusuf, “Extended block hybrid backward differentiation formula for second order fuzzy differen- tial equations using legendre polynomial as basis function”, J. Sci. Tech. Math. Educ. 16 (2020) 100. [17] S. Al-Refai, M. I. Syam, & M. Al-Refai, “Optimization of one-step block method for solving second-order fuzzy initial value problems”, Complex. 2021 (2021) 1. https://doi.org/10.1155/2021/6650413 [18] B. Bede, “Fuzzy sets”, in Mathematics of Fuzzy Sets and Fuzzy Logic, Springer (2013) 1. https://doi.org/10.1007/978-3-642-35221-8 1 9 Hussain et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1087 10 [19] S. S. Devi & K. Ganesan, “An approximate solution by fuzzy Taylor’s method”, Int. J. Pure Appl. Math. 113 (2017) 236. [20] J. D. Lambert, Computational methods in ordinary differen- tial equations, New York: John Wiley and Sons Inc, (1973). https://doi.org/10.1002/zamm.19740540726 [21] J. C. Butcher, Numerical methods for ordinary differential equations, New York: John Wiley and Sons Ltd (2008). https://doi.org/10.1002/9781119121534 [22] E. Süli, Numerical solution of ordinary differential equations, Mathemat- ical Institute, University of Oxford (2010). [23] J. O. Ehigie, S. A. Okunuga & A. B. Sofoluwe, “3-point block methods for direct integration of general second-order ordinary differential equa- tions”, Advances in Numerical Analysis (2011). [24] M. O. Ogunniran, “A class of block multi-derivative numerical techniques for singular advection equations”, J. Nig. Soc. Phy. Sci. 1 (2019) 62. https://doi.org/10.46481/jnsps.2019.12 [25] R. Abdulganiy, O. Akinfenwa, O. Yusuff, O. Enobabor, & S. Okunuga, “Block third derivative trigonometrically-fitted methods for stiff and periodic problems”, J. Nig. Soc. Phy. Sci. 2 (2020) 12. https://doi.org/10.46481/jnsps.2020.33 [26] Z. Omar & O. Adeyeye, “Numerical solution of first order initial value problems using a self-starting implicit two-step Obrechkoff-type block method”, J. Math. Stat. 12 (2016) 127. http://doi.org/10.3844/jmssp.2016.127.134 [27] Z. Zlatev, K. Georgiev, & I. Dimov, “Studying absolute stability properties of the Richardson Extrapolation combined with explicit Runge–Kutta methods”, Comput. Math. with Appl. 67 (2014) 2294. https://doi.org/10.1016/j.camwa.2014.02.025 [28] A. O. Adesanya, D. M. Udoh, & A. M. Ajileye, “A new hybrid block method for the solution of general third order initial value problems of ordinary differential equations”, Int. J. Pure Appl. Math. 86 (2013) 365. http://dx.doi.org/10.12732/ijpam.v86i2.11 [29] J. O. Kuboye & Z. Omar, “Numerical solution of third order ordinary dif- ferential equations using a seven-step block method”, Int. J. Math. Anal. 9 (2014) 743. http://dx.doi.org/10.12988/ijma.2015.5125 [30] A. F. Jameel, A. Saaban, & H. H. Zureigat, “Numerical solution of second-order fuzzy nonlinear two-point boundary value problems using combination of finite difference and Newton’s methods”, Neural. Com- put. Appl. 30 (2018) 3167. https://doi.org/10.1007/s00521-017-2893-z 10