EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 16, No. 2, 2023, 1236-1259 ISSN 1307-5543 – ejpam.com Published by New York Business Global Solving Fuzzy System of Boundary Value Problems by Homotopy Perturbation Method with Green’s Function Waleed Al-Hayani1,∗, Mahasin Thabet Younis1 1 Department of Mathematics, College of Computer Science and Mathematics, University of Mosul, Iraq. Abstract. In this research, we successfully demonstrated the use of the homotopy perturbation method with Green’s function to find approximate solutions for the fuzzy system of boundary value problems. Our results showcase the effectiveness of this method in providing accurate and reliable solutions. The consistent way to reduce the size of the computation gives to reach the exact solution is one of the best methods adopted to determine the behavior of the solution directly in order to determine the approximate solution analytically, Finally, the problems that have been addressed confirmed the validity of the method applied analytically in this research using comparison with some numerical problems. 2020 Mathematics Subject Classifications: 03E72, 34A07, 34B05, 34B15 Key Words and Phrases: Fuzzy system, Boundary value problems, Homotopy perturbation method, He’s polynomials, Green’s function 1. Introduction A multitude of dynamic processes can be represented mathematically as systems of partial or ordinary differential equations. Differential equations have been widely rec- ognized as a successful modeling approach. Meanwhile, fuzzy group theory serves as a potent tool in handling unknown variables and processing fuzzy or subjective information in mathematical models [1, 5, 9, 18]. It may be necessary to use such mathematical mod- eling, and to use fuzzy differential equations (FDEs) that involve a system of elementary value problems. Fuzzy systems of differential equations appear when modeling these prob- lems are incomplete and questionable in the case of inaccurate data. Fuzzy differential equation systems are applicable mathematical models [6]. Dynamic models of dynamic systems with uncertainties or ambiguous models, represented by fuzzy differential equa- tions, are widely utilized across a spectrum of fields, such as geometry and population modeling. These models have proven to be versatile and applicable in a diverse range of ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v16i2.4764 Email addresses: waleedalhayani@uomosul.edu.iq (Waleed Al-Hayani), mahasin thabet@uomosul.edu.iq (Mahasin Younis) https://www.ejpam.com 1236 © 2023 EJPAM All rights reserved. W. Al-Hayani, M T. Younis / Eur. J. Pure Appl. Math, 16 (2) (2023), 1236-1259 1237 applications, making them a valuable tool for data analysis and modeling [24]. In recent years, the homotopy perturbation method (HPM) has been studied and success- fully applied by many engineers, scientists and researchers, and used to solve the integral equations and differential equations. This HPM was introduced first by Dr. Ji Huan He in 1998 [15, 16]. The HPM is a coupling of the homotopy method, a basic concept of topology, and the classical perturbation technique. This coupling will provide with a suitable way to obtain approximate or analytic solution for different problems arising in various scientific fields [14]. The HPM boasts several advantages, including the ability to obtain accurate approximate solutions with only a few iterations and a rapid convergence of the solution series in many instances. This method has proven to be a valuable tool in solving a wide range of problems, particularly in the realm of integral equations [20]. The integral equations can be used to describe a wide variety of problems in science and engineering [3]. The HPM is a sophisticated and effective approach for solving system of linear and nonlinear differential and integral equations [2]. A system of linear and non- linear differential equations may also be solved using this approach due to the complexity of finding exact solutions for nonlinear differential equations, any perturbative method that meets certain criteria is deemed acceptable. The HPM offers an analytical solution by utilizing the initial conditions and is noteworthy for its ability to produce highly accu- rate approximate solutions with only a small number of terms required [4]. In most cases, only approximate solutions, either numerical or analytical, can be obtained for nonlinear problems. The numerous branches of science have a vast array of nonlinear problems for which exact solutions cannot be found. As a result, numerous analytical and numerical approximations have been studied and developed to tackle these complex equations. Thus, alternative methods become crucial for solving these nonlinear equations effectively [14]. The concept of fuzziness has experienced a surge in popularity and application across various fields since its inception, including learning theory, automata, decision-making processes, algorithms, pattern classification, and linguistics. The idea behind fuzzy con- cepts involves converting differential equations into differential inclusions and accepting the solution as the α−cut of the fuzzy solution, leading to the transformation of the fuzzy differential equation (FDE) into a equivalent fuzzy integral equation [6]. In the last years, the researchers Visuvasam et al. [23], Shanthi et al. [19], Vijayalakshmi et al. [21] and Vijayalakshmi et al. [22] have applied HPM to solve the mathematical models as (porous rotating disk electrodes, ECE reactions at rotating disk electrodes, magnetohydrodynamic, etc.). The main objective of the fuzzy HPM with Green’s function, is to develop a more effec- tive and efficient method for solving linear and non-linear fuzzy system of boundary value problems (FSBVPs) 2∑ i=0 ai (x) ỹ (i) 1 (x, α) + 1∑ i=0 bi (x) ỹ (i) 2 (x, α) +N1 [ỹ1 (x, α) , ỹ2 (x, α)] = f̃1 (x, α) , 2∑ i=0 ci (x) ỹ (i) 2 (x, α) + 1∑ i=0 di (x) ỹ (i) 1 (x, α) +N2 [ỹ1 (x, α) , ỹ2 (x, α)] = f̃2 (x, α) , x ∈ [0, 1] (1) W. Al-Hayani, M T. Younis / Eur. J. Pure Appl. Math, 16 (2) (2023), 1236-1259 1238 with the boundary conditions ỹ1 (0, α) = ỹ1 (1, α) = (0, 0) , ỹ2 (0, α) = ỹ2 (1, α) = (0, 0) , (2) where ỹ1 (x, α) , ỹ2 (x, α) , f̃1 (x, α) , f̃2 (x, α) and ai (x) , bi (x) , ci (x) , di (x) for i = 0, 1, 2 are analytic functions, and we are seeking the solutions ỹ1 (x, α) , ỹ2 (x, α) satisfying Equa- tion 1. We assume that for every f̃1 (x, α) , f̃2 (x, α) Equation 1 has one and only one solution. By combining the advantages of the HPM and Green’s function, this study aims to provide a new approach that leverages the strengths of both methods to obtain accurate approximate solutions for fuzzy ODEs. This research can help to deepen our understand- ing of fuzzy system of ODEs and provide valuable insights into the development of new techniques for solving these complex problems. 2. The fundamental concept of HPM Now, we demonstrate the fundamental concept of the HPM [8, 13, 14], we assume the following non-linear differential equation. D (y)− g (r) = 0, r ∈ Ω (3) with the condition B ( y, ∂y ∂n ) = 0, r ∈ Γ (4) where D is a general differential operator, which can be divided into linear parts L and N nonlinear term of Equation 4, B is a boundary operator, g(r) is a known analytic function, Γ is the boundary of the domain Ω. can be rewritten as L(y) +N(y)− g(r) = 0. (5) By the homotopy technique, we construct a homotopy v(r, p) : Ω × [0, 1] → R which satisfies H(v, p) = (−p)[L(v)− L(v0)] + p[A(v)− g(r)] = 0, p ∈ [0, 1], r ∈ Ω) (6) or H(v, p) = L(v)− L(v0) + pL(v0) + p[N(v)− g(r)] = 0, (7) where p ∈ [0, 1] is an embedding parameter, v0 is an initial approximation of Equation 6 which satisfies the boundary conditions. Obviously, from Equation 7 we have H(v, 0) = L(v)− L(v0) = 0, (8) H(v, 1) = D(v)− g(r) = 0. (9) W. Al-Hayani, M T. Younis / Eur. J. Pure Appl. Math, 16 (2) (2023), 1236-1259 1239 The changing process of p from zero to unity is just that of y(r, p) from y0(r) to y(r). In topology, this is called deformation, and L(y)− L(v0), B(y)− f(r) are called homptopic. In this paper, the authors will first use the imbedding parameter p as a “small pa- rameter”, and assume that the solution of Equation 7 can be written as a power series in p : y = ∞∑ n=0 pnyn = y0 + py1 + p2y2 + · · · , (10) and the nonlinear term Ny can be decomposed as Ny = ∞∑ n=0 pnHn (y) , (11) where the Hn are He’s polynomials of y0, y1, . . . , yn and are calculated by the definitional formula [11] Hn (y0, y1, . . . , yn) = 1 n! ∂n ∂pn [ N ( ∞∑ i=0 piyi )] p=0 , n = 0, 1, 2, . . . . Setting p = 1 results in the approximate solution of Equation 3 y = lim p→1 ∞∑ n=0 pnyn = y0 + y1 + y2 + y3 + · · · . (12) The convergence of the series solution (10) is given in [12]. 3. HPM with Green’s function The Green’s function G(x, ξ) corresponding to linear operator d2 dx2 ỹi (x, α) , i = 1, 2 in Eq. (1) is G(x, ξ) =  C1 (ξ)x+ C2 (ξ) , 0 ≤ ξ ≤ x, C3 (ξ)x+ C4 (ξ) , x ≤ ξ ≤ 1, (13) where x ̸= ξ, C1 (ξ) , C2 (ξ) , C3 (ξ) and C4 (ξ) are linearly independent solutions of d2 dx2 ỹi (x, α) with a2 (x) = 1 and c2 (x) = 1. We follow three conditions to find the Green’s function of the second-order BVP. – G(x, ξ) is continuous at x = ξ : [C1 (ξ)x+ C2 (ξ)]x=ξ = [C3 (ξ)x+ C4 (ξ)]x=ξ (14) – ∂ ∂x G(x, ξ) is discontinuous at x = ξ :{ ∂ ∂x [C1 (ξ)x+ C2 (ξ)]− ∂ ∂x [C3 (ξ)x+ C4 (ξ)] } x=ξ = −1 (15) W. Al-Hayani, M T. Younis / Eur. J. Pure Appl. Math, 16 (2) (2023), 1236-1259 1240 – G(x, ξ) satisfies the homogeneous boundary conditions (2): [C1 (ξ)x+ C2 (ξ)]x=0 = 0, [C3 (ξ)x+ C4 (ξ)]x=1 = 0. (16) By solving the above three equations (14)-(16) to find C1 (ξ) , C2 (ξ) , C3 (ξ) and C4 (ξ), therefore, the Green’s function (13) of Eqs. (1) and (2) can be written as G(x, ξ) =  ξ (1− x) , 0 ≤ ξ ≤ x, x (1− ξ) , x ≤ ξ ≤ 1. (17) Now, we costruct the HPM with Green’s function, we have ỹ1 (x, α) = h̃1 (x, α)− ∫ 1 0 G (x, ξ) f̃1 (ξ, α) dξ +p ∫ 1 0 G (x, ξ) 1∑ i=0 ai (ξ) ỹ (i) 1 (ξ, α) dξ +p ∫ 1 0 G (x, ξ) 1∑ i=0 bi (ξ) ỹ (i) 2 (ξ, α) dξ +p ∫ 1 0 G (x, ξ)N1 [ỹ1 (ξ, α) , ỹ2 (ξ, α)] dξ, ỹ2 (x, α) = h̃2 (x, α)− ∫ 1 0 G (x, ξ) f̃2 (ξ, α) dξ +p ∫ 1 0 G (x, ξ) 1∑ i=0 ci (ξ) ỹ (i) 2 (ξ, α) dξ, +p ∫ 1 0 G (x, ξ) 1∑ i=0 di (ξ) ỹ (i) 1 (ξ, α) dξ +p ∫ 1 0 G (x, ξ)N2 [ỹ1 (ξ, α) , ỹ2 (ξ, α)] dξ, (18) where h̃i (x, α) is the solution of y′′i (x, α) = 0, i = 1, 2 with the boundary conditions (2). W. Al-Hayani, M T. Younis / Eur. J. Pure Appl. Math, 16 (2) (2023), 1236-1259 1241 Substituting Eqs. 10 and 11 into Equation 18, we obtain ∞∑ n=0 pnỹ1,n (x, α) = h̃1 (x, α)− ∫ 1 0 G (x, ξ) f̃1 (ξ, α) dξ +p ∫ 1 0 G (x, ξ) 1∑ i=0 ai (ξ) ∞∑ n=0 pnỹ (i) 1,n (ξ, α) dξ +p ∫ 1 0 G (x, ξ) 1∑ i=0 bi (ξ) ∞∑ n=0 pnỹ (i) 2,n (ξ, α) dξ +p ∫ 1 0 G (x, ξ) ∞∑ n=0 pnH̃1,ndξ ỹ2 (x, α) = h̃2 (x, α)− ∫ 1 0 G (x, ξ) f̃2 (ξ, α) dξ +p ∫ 1 0 G (x, ξ) 1∑ i=0 ci (ξ) ∞∑ n=0 pnỹ (i) 2,n (ξ, α) dξ, +p ∫ 1 0 G (x, ξ) 1∑ i=0 di (ξ) ∞∑ n=0 pnỹ (i) 1,n (ξ, α) dξ +p ∫ 1 0 G (x, ξ) ∞∑ n=0 pnH̃2,ndξ. (19) Using the fuzzy HPM, according to Equation 19, the lower iterations (L) are then deter- mined in the following recursive way: p0 :  y 1,0 (x, α) = h1 (x, α)− ∫ 1 0 G (x, ξ) f 1 (ξ, α) dξ, y 2,0 (x, α) = h2 (x, α)− ∫ 1 0 G (x, ξ) f 2 (ξ, α) dξ, pn :  y 1,n+1 (x, α) = ∫ 1 0 G (x, ξ) 1∑ i=0 ai (ξ) y (i) 1,n (ξ, α) dξ, + ∫ 1 0 G (x, ξ) [ 1∑ i=0 bi (ξ) y (i) 2,n (ξ, α) +H1,n ] dξ, y 2,n+1 (x, α) = ∫ 1 0 G (x, ξ) 1∑ i=0 ci (ξ) y (i) 2,n (ξ, α) dξ, + ∫ 1 0 G (x, ξ) [ 1∑ i=0 di (ξ) y (i) 1,n (ξ, α) +H2,n ] dξ, , n ≥ 1 (20) W. Al-Hayani, M T. Younis / Eur. J. Pure Appl. Math, 16 (2) (2023), 1236-1259 1242 and the upper iterations (U) are p0 : { y1,0 (x, α) = h1 (x, α)− ∫ 1 0 G (x, ξ) f1 (ξ, α) dξ, y2,0 (x, α) = h2 (x, α)− ∫ 1 0 G (x, ξ) f2 (ξ, α) dξ, pn :  y1,n+1 (x, α) = ∫ 1 0 G (x, ξ) 1∑ i=0 ai (ξ) y (i) 1,n (ξ, α) dξ, + ∫ 1 0 G (x, ξ) [ 1∑ i=0 bi (ξ) y (i) 2,n (ξ, α) +H1,n ] dξ, y2,n+1 (x, α) = ∫ 1 0 G (x, ξ) 1∑ i=0 ci (ξ) y (i) 2,n (ξ, α) dξ, + ∫ 1 0 G (x, ξ) [ 1∑ i=0 di (ξ) y (i) 1,n (ξ, α) +H2,n ] dξ. , n ≥ 1 (21) By solving the system (20) and (21), we obtain the iterations ỹi,0 (x, α), ỹi,1 (x, α) , . . . , ỹi,n (x, α) , i = 1, 2. Thus, the approximated solution in a series form is given by ỹi (x, α) = lim p→1 ∞∑ n=0 pnỹi,n (x, α) = ∞∑ n=0 ỹi,n (x, α) . 4. Basic Concepts Definition 1. [10] A fuzzy number ũ : R → [0, 1] satisfying the properties: (i) ũ is normal, if ũ (x0) = 1, x0 ∈ R. (ii) ũ is fuzzy convex set if ũ (λx+ (1− λ) y) ≥ min {ũ (x) , ũ (y)}, ∀ x, y ∈ R, λ ∈ [0, 1] . (iii) ũ is upper semi-continuous on R. (iv) ũ (x),closure of ũ (x) x ∈ R and is called a compact set. The set of all fuzzy numbers is denoted by RF . For 0 < α ≤ 1, denote [ũ]α = {x ∈ R, ũ (x) ≥ α} and [ũ]0 = {x ∈ R, ũ (x) > 0}. Then, it is well-known that for any α ∈ [0, 1] , [ũ]α is a bounded closed interval. A triangular fuzzy numbers are very popular and denoted by A = (a1, a2, a3) and defined by µA (x) =  0, x < a1 x− a1 a2−a1 , a1 ≤ x ≤ a2 a3 − x a3 − a2 , a2 ≤ x ≤ a3 0, x > a3 W. Al-Hayani, M T. Younis / Eur. J. Pure Appl. Math, 16 (2) (2023), 1236-1259 1243 where a1, a3 > 0. Definition 2. [10, 24] An arbitrary fuzzy number ũ in the parametric form is represented by an ordered pair of functions ( y (x, α) , y (x, α) ) , 0 ≤ α ≤ 1; which satisfy the following requirements . (i) y (x, α) is a bounded left continuous non-decreasing function over [0, 1] . (ii) y (x, α)is a bounded left continuous non-increasing function over [0, 1] . (iii) y (x, α) ≤ y (x, α) , 0 ≤ α ≤ 1. Definition 3. [24] For arbitrary fuzzy numbers ũ (x, α) = (u (x, α) , u (x, α)) and ṽ (x, α) = (v (x, α) , v (x, α)) , we define the distance between ũ (x, α) and ṽ (x, α) by the quantity D (ũ (x, α) , ṽ (x, α)) = sup α∈[0,1] max {|u (x, α)− v (x, α)| , |u (x, α)− v (x, α)|} 5. Applications In this section, we demonstrate the exceptional precision and accuracy of the results obtained through this method. By showcasing the solutions of two distinct problems.we give a comparison the approximate solution by the (HPM with Green’s function) with the exact solution given. Moreover, the absolute errors between them are defined as follows: AEL = ∣∣∣y i,E (x, α)− y i,n (x, α) ∣∣∣ AEU = ∣∣yi,E (x, α)− yi,n (x, α) ∣∣ , i = 1, 2 n = 1, 2, . . . where y i,E (x, α) : Lower exact solution, y i,n (x, α) : Lower approximate solution (HPM), yi,E (x, α) : Upper exact solution, yi,n (x, α) : Upper approximate solution (HPM). The computations associated with the problems were performed using a Maple 18 package with a precision of 20 digits. Problem 1. We first consider the linear FSBVPs [7, 9, 17]:{ ỹ′′1 (x, α) + (2x− 1) ỹ′1 (x, α) + cos (πx) ỹ′2 (x, α) = f̃1 (x, α) , ỹ′′2 (x, α) + xỹ1 (x, α) = f̃2 (x, α) , (22) W. Al-Hayani, M T. Younis / Eur. J. Pure Appl. Math, 16 (2) (2023), 1236-1259 1244 with the boundary conditions ỹ1 (0, α) = ỹ1 (1, α) = (0, 0) , ỹ2 (0, α) = ỹ2 (1, α) = (0, 0) , (23) where f̃1 (x, α) = [ f 1 (x, α) , f1 (x, α) ] and f̃2 (x, α) = [ f 2 (x, α) , f2 (x, α) ] are given by f 1 (x, α) = −α sin (πx)π2 + [(2απ + α+ 1)x] cos (πx) , − [ απ + 1 2 α+ 1 2 ] cos (πx) , f1 (x, α) = (α− 2)π2 sin (πx) + [(4π − 2απ − 4α+ 6)x] cos (πx) , + [απ + 2α− 2π − 3] cos (πx) , f 2 (x, α) = xα sin (πx) + α+ 1, f2 (x, α) = (2− α)x sin (πx)− 4α+ 6. (24) The exact solutions of the FSBVPs (22) and (23) are ỹ1,E (x, α) = [α sin (πx) , (2− α) sin (πx)] , ỹ2,E (x, α) = [( α 2 + 1 2 )( x2 − x ) , (3− 2α) ( x2 − x )] . Constructing the HPM with Green’s function, we have ỹ1 (x, α) = h̃1 (x, α)− ∫ 1 0 G (x, ξ) f̃1 (ξ, α) dξ +p ∫ 1 0 G (x, ξ) [(2ξ − 1) ỹ′1 (ξ, α) + cos (πx) ỹ′2 (ξ, α)] dξ, ỹ2 (x, α) = h̃2 (x, α)− ∫ 1 0 G (x, ξ) f̃2 (ξ, α) dξ +p ∫ 1 0 G (x, ξ) [ξỹ1 (ξ, α)] dξ, (25) where h̃i (x, α) is the solution of y′′i (x, α) = 0, i = 1, 2 with the boundary conditions (23) and G(x, ξ) is the Green’s function given by (17). Substituting Equation 10 into Equation 25, we obtain ∞∑ n=0 pnỹ1,n (x, α) = h̃1 (x, α)− ∫ 1 0 G (x, ξ) f̃1 (ξ, α) dξ +p ∫ 1 0 G (x, ξ) [ (2ξ − 1) ∞∑ n=0 pnỹ′1,n (ξ, α) + cos (πx) ∞∑ n=0 pnỹ′2,n (ξ, α) ] dξ, ∞∑ n=0 pnỹ2,n (x, α) = h̃2 (x, α)− ∫ 1 0 G (x, ξ) f̃2 (ξ, α) dξ +p ∫ 1 0 G (x, ξ) [ ξ ∞∑ n=0 pnỹ1,n (ξ, α) ] dξ. (26) W. Al-Hayani, M T. Younis / Eur. J. Pure Appl. Math, 16 (2) (2023), 1236-1259 1245 Using the fuzzy HPM, according to Equation 26, the lower iterations (L) are then deter- mined in the following recursive way: p0 :  y 1,0 (x, α) = h1 (x, α)− ∫ 1 0 G (x, ξ) f 1 (ξ, α) dξ, y 2,0 (x, α) = h2 (x, α)− ∫ 1 0 G (x, ξ) f 2 (ξ, α) dξ, pn :  y 1,n+1 (x, α) = ∫ 1 0 G (x, ξ) [ (2ξ − 1) y′ 1,n (ξ, α) + cos (πx) y′ 2,n (ξ, α) ] dξ, y 2,n+1 (x, α) = ∫ 1 0 G (x, ξ) [ ξy 1,n (ξ, α) ] dξ, (27) and the upper iterations (U) are p0 : { y1,0 (x, α) = h1 (x, α)− ∫ 1 0 G (x, ξ) f1 (ξ, α) dξ, y2,0 (x, α) = h2 (x, α)− ∫ 1 0 G (x, ξ) f2 (ξ, α) dξ, pn : { y1,n+1 (x, α) = ∫ 1 0 G (x, ξ) [ (2ξ − 1) y′1,n (ξ, α) + cos (πx) y′2,n (ξ, α) ] dξ, y2,n+1 (x, α) = ∫ 1 0 G (x, ξ) [ ξy1,n (ξ, α) ] dξ. (28) By solving the system (27) and (28), we obtain the iterations ỹi,0 (x, α), ỹi,1 (x, α) , . . . , ỹi,n (x, α) , i = 1, 2. Thus, the approximated solution in a series form is given by ỹi (x, α) = 5∑ n=0 ỹi,n (x, α) . We compare the approximated solution using the HPM (n = 5) with the exact solution and the error produced in Tables 1–4. Table 5 demonstrate the maximum errors on the interval 0 ≤ x ≤ 1. W. Al-Hayani, M T. Younis / Eur. J. Pure Appl. Math, 16 (2) (2023), 1236-1259 1246 Table 1: Comparison of the numerical results for y 1 (x, α) (Problem 1) α x y 1,E (x, α) y 1,5 (x, α) AEL 0.25 0.1 0.0772542 0.0772587 4.531E − 06 0.3 0.2022542 0.2022542 6.846E − 06 0.5 0.2500000 0.2500069 6.923E − 06 0.7 0.2022542 0.2022610 6.771E − 06 0.9 0.0772542 0.0772589 4.728E − 06 0.50 0.1 0.1545084 0.1545163 7.849E − 06 0.3 0.4045084 0.4045203 1.185E − 05 0.5 0.5000000 0.5000119 1.198E − 05 0.7 0.4045084 0.4045202 1.172E − 05 0.9 0.1545084 0.1545166 8.188E − 06 0.75 0.1 0.2317627 0.2317739 1.116E − 05 0.3 0.6067627 0.6067796 1.685E − 05 0.5 0.7500000 0.7500170 1.704E − 05 0.7 0.6067627 0.6067794 1.667E − 05 0.9 0.2317627 0.2317743 1.164E − 05 1.00 0.1 0.3090169 0.3090314 1.448E − 05 0.3 0.8090169 0.8090388 2.186E − 05 0.5 1.0000000 1.0000221 2.210E − 05 0.7 0.8090169 0.8090386 2.162E − 05 0.9 0.3090169 0.3090321 1.510E − 05 Table 2: Comparison of the numerical results for y1 (x, α) (Problem 1) α x y1,E (x, α) y1,5 (x, α) AEU 0.25 0.1 0.5407797 0.5408069 2.716E − 05 0.3 1.4157797 1.4158207 4.101E − 05 0.5 1.75000000 1.7500414 4.147E − 05 0.7 1.4157797 1.4158203 4.056E − 05 0.9 0.5407797 0.5408080 2.834E − 05 0.50 0.1 0.4635254 0.4635484 2.294E − 05 0.3 1.2135254 1.2135601 3.463E − 05 0.5 1.5000000 1.5000350 3.501E − 05 0.7 1.2135254 1.2135597 3.425E − 05 0.9 0.4635254 0.4635494 2.393E − 05 0.75 0.1 0.3862712 0.3862899 1.871E − 05 0.3 1.0112712 1.0112994 2.824E − 05 0.5 1.2500000 1.2500285 2.855E − 05 0.7 1.0112712 1.0112991 2.793E − 05 0.9 0.3862712 0.3862907 1.952E − 05 1.00 0.1 0.3090169 0.3090314 1448E − 05 0.3 0.8090169 0.8090388 2.186E − 05 0.5 1.0000000 1.0000221 2.210E − 05 0.7 0.8090169 0.8090386 2.162E − 05 0.9 0.3090169 0.3090321 1.510E − 05 W. Al-Hayani, M T. Younis / Eur. J. Pure Appl. Math, 16 (2) (2023), 1236-1259 1247 Table 3: Comparison of the numerical results for y 2 (x, α) (Problem 1) α x y 2,E (x, α) y 2,5 (x, α) AEL 0.25 0.1 −0.056250 −0.0562509 9.101E − 07 0.3 −0.131250 −0.1312525 2.542E − 06 0.5 −0.156250 −0.1562534 3.486E − 06 0.7 −0.131250 −0.1312532 3.256E − 06 0.9 −0.056250 −0.0562514 1.452E − 06 0.50 0.1 −0.067500 −0.6750158 1.583E − 05 0.3 −0.157500 −0.1575044 4.423E − 06 0.5 −0.187500 −0.1875060 6.065E − 06 0.7 −0.157500 −.01575056 5.665E − 06 0.9 −0.067500 −0.0675025 2.527E − 06 0.75 0.1 −0.787500 −0.7875225 2.256E − 06 0.3 −0.183750 −0.1837563 6.304E − 06 0.5 −0.218750 −0.2187586 8.644E − 06 0.7 −0.183750 −0.1837580 8.074E − 06 0.9 −0.787500 −0.7875360 3.601E − 06 1.00 0.1 −0.900000 −0.9000292 2.929E − 06 0.3 −0.210000 −0.2100081 8.186E − 06 0.5 −0.250000 −0.2500112 1.122E − 05 0.7 −0.210000 −0.2100104 1.048E − 05 0.9 −0.900000 −0.9000467 4.675E − 06 Table 4: Comparison of the numerical results for y2 (x, α) (Problem 1) α x y2,E (x, α) y2,5 (x, α) AEU 0.25 0.1 −0.2250 −0.2250054 5.048E − 06 0.3 −0.5250 −0.5250153 1.531E − 05 0.5 −0.6250 −0.6250210 2.100E − 05 0.7 −0.5250 −0.5250196 1.961E − 05 0.9 −0.2250 −0.2250087 8.750E − 06 0.50 0.1 −0.1800 −0.1800046 4.631E − 06 0.3 −0.3200 −0.3200090 9.056E − 06 0.5 −0.5000 −0.5000177 1.774E − 05 0.7 −0.4200 −.04200165 1.657E − 05 0.9 −0.1800 −0.1800073 7.392E − 06 0.75 0.1 −0.1350 −0.1350037 3.780E − 06 0.3 −0.3150 −0.3150105 1.056E − 05 0.5 −0.3750 −0.3750144 1.448E − 05 0.7 −0.3150 −0.3150135 1.352E − 05 0.9 −0.1350 −0.1350060 6.034E − 06 1.00 0.1 −0.9000 −0.9000292 2.929E − 06 0.3 −0.2100 −0.2100081 8.186E − 06 0.5 −0.2500 −0.2500112 1.122E − 05 0.7 −0.2100 −0.2100104 1.048E − 05 0.9 −0.9000 −0.9000467 4.675E − 06 W. Al-Hayani, M T. Younis / Eur. J. Pure Appl. Math, 16 (2) (2023), 1236-1259 1248 Table 5: Maximum errors (Problem 1) ỹ1 (x, α) ỹ2 (x, α) α n y 1 (x, α) y1 (x, α) y 2 (x, α) y2 (x, α) 0.25 5 6.957E − 06 4.167E − 05 3.546E − 06 2.136E − 05 0.50 5 1.204E − 05 3.519E − 05 6.169E − 06 1.804E − 05 0.75 5 1.712E − 05 2.870E − 05 8.793E − 06 1.473E − 05 1.00 5 2.221E − 05 2.221E − 05 1.141E − 05 1.141E − 05 In Figures 1–4, a very good agreement is shown between the exact solution ( ỹi,E(x, 0.5), i = 1, 2 ) with a continuous line and the approximate solution by the HPM with Green’s function ( ỹi,5(x, 0.5), i = 1, 2 ) with the symbol o. We present the contour plot in 2D on the (x, α)-plane for the exact solutions ( ỹi,E(x, 0.5), i = 1, 2 ) with a continuous line and the approximate solution by the HPM with Green’s function (ỹi,5(x, 0.5), i = 1, 2) with the symbol o in Figures 5–8. Figure 1: Plot of y 1,E (x, 0.5) and y 1,5 (x, 0.5). Figure 2: Plot of y1,E(x, 0.5) and y1,5(x, 0.5). W. Al-Hayani, M T. Younis / Eur. J. Pure Appl. Math, 16 (2) (2023), 1236-1259 1249 Figure 3: Plot of y 2,E (x, 0.5) and y 2,5 (x, 0.5). Figure 4: Plot of y2,E(x, 0.5) and y2,5(x, 0.5). Figure 5: The contour plot in 2D (x, α)-plane for y 1,E (x, α) and y 1,5 (x, α). W. Al-Hayani, M T. Younis / Eur. J. Pure Appl. Math, 16 (2) (2023), 1236-1259 1250 Figure 6: The contour plot in 2D (x, α)-plane for y1,E(x, α) and y1,5(x, α). Figure 7: The contour plot in 2D (x, α)-plane for y 2,E (x, α) and y 2,5 (x, α). Figure 8: The contour plot in 2D (x, α)-plane for y2,E(x, α) and y2,5(x, α). W. Al-Hayani, M T. Younis / Eur. J. Pure Appl. Math, 16 (2) (2023), 1236-1259 1251 Problem 2 Now, we turn to a non-linear FSBVPs [7, 9, 17]:{ ỹ′′1 (x, α) + xỹ′1 (x, α) + cos (πx) ỹ′2 (x, α) = f̃1 (x, α) , ỹ′′2 (x, α) + xỹ′1 (x, α) + xy21 (x, α) = f̃2 (x, α) , (29) with the boundary conditions ỹ1 (0, α) = ỹ1 (1, α) = (0, 0) , ỹ2 (0, α) = ỹ2 (1, α) = (0, 0) , (30) where f̃1 (x, α) = [ f 1 (x, α) , f1 (x, α) ] and f̃2 (x, α) = [ f 2 (x, α) , f2 (x, α) ] are given by f 1 (x, α) = α2 sin (x) + (α2x2 − α2x+ 2α2) cos (x) +(−α3x− x+ α3 2 + 1 2 ) cos (πx) , f1 (x, α) = (3− 2α) sin (x) + (−2αx2 + 3x2 + 2αx− 3x− 4α+ 6) cos (x) +(2α3x− 4x− α3 + 2) cos (πx) , f 2 (x, α) = ( α4x3 − 2α4x2 + α4x ) sin (x)2 + α2x sin(x) +(α2x2 − α2x) cos (x)− α3 − 1, f2 (x, α) = + ( 4α2 − 12α+ 9 ) x3 sin (x)2 + (−8α2 + 24α− 18)x2 sin (x)2 + ( 4α2 − 12α+ 9 ) x sin (x)2 + (−2α+ 3)x sin (x) +(−2α+ 3)x2 cos (x) + (2α− 3)x cos (x)− 4 + 2α3. (31) The exact solutions of the FSBVPs (29) and (30) are ỹ1,E (x, α) = [ α2(x− 1) sin (x) , (3− 2α) (x− 1) sin (x) ] , ỹ2,E (x, α) = [ (α3 + 1) 2 ( x− x2 ) , ( 2− α3 ) ( x− x2 )] . Constructing the HPM with Green’s function, we have ỹ1 (x, α) = h̃1 (x, α)− ∫ 1 0 G (x, ξ) f̃1 (ξ, α) dξ +p ∫ 1 0 G (x, ξ) [ξỹ′1 (ξ, α) + cos (πx) ỹ′2 (ξ, α)] dξ, ỹ2 (x, α) = h̃2 (x, α)− ∫ 1 0 G (x, ξ) f̃2 (ξ, α) dξ +p ∫ 1 0 G (x, ξ) [ ξỹ′1 (ξ, α) + ξỹ21 (ξ, α) ] dξ, (32) where h̃i (x, α) is the solution of y′′i (x, α) = 0, i = 1, 2 with the boundary conditions (30) and G(x, ξ) is the Green’s function given by (17). Substituting Eqs. (10) and (11) into W. Al-Hayani, M T. Younis / Eur. J. Pure Appl. Math, 16 (2) (2023), 1236-1259 1252 Equation 32, we obtain ∞∑ n=0 pnỹ1,n (x, α) = h̃1 (x, α)− ∫ 1 0 G (x, ξ) f̃1 (ξ, α) dξ +p ∫ 1 0 G (x, ξ) [ ξ ∞∑ n=0 pnỹ′1,n (ξ, α) + cos (πx) ∞∑ n=0 pnỹ′2,n (ξ, α) ] dξ, ∞∑ n=0 pnỹ2,n (x, α) = h̃2 (x, α)− ∫ 1 0 G (x, ξ) f̃2 (ξ, α) dξ +p ∫ 1 0 G (x, ξ) [ ξ ∞∑ n=0 pnỹ′1,n (ξ, α) + ξ ∞∑ n=0 pnH1,n ] dξ. (33) Using the fuzzy HPM, according to Equation 33, the lower iterations (L) are then deter- mined in the following recursive way: p0 :  y 1,0 (x, α) = h1 (x, α)− ∫ 1 0 G (x, ξ) f 1 (ξ, α) dξ, y 2,0 (x, α) = h2 (x, α)− ∫ 1 0 G (x, ξ) f 2 (ξ, α) dξ, pn :  y 1,n+1 (x, α) = ∫ 01G (x, ξ) [ ξy′ 1,n (ξ, α) + cos (πx) y′ 2,n (ξ, α) ] dξ, y 2,n+1 (x, α) = ∫ 01G (x, ξ) [ ξy′ 1,n (ξ, α) + ξH1,n ] dξ, (34) and the upper iterations (U) are p0 : { y1,0 (x, α) = h1 (x, α)− ∫ 1 0 G (x, ξ) f1 (ξ, α) dξ, y2,0 (x, α) = h2 (x, α)− ∫ 1 0 G (x, ξ) f2 (ξ, α) dξ, pn : { y1,n+1 (x, α) = ∫ 1 0 G (x, ξ) [ ξy′1,n (ξ, α) + cos (πx) y′2,n (ξ, α) ] dξ, y2,n+1 (x, α) = ∫ 1 0 G (x, ξ) [ ξy′1,n (ξ, α) + ξH1,n ] dξ. (35) where the nonlinear terms defined by the series y2 1,n (ξ, α) = ∑∞ n=0H1,n, y21,n (ξ, α) = ∑∞ n=0H1,n, and the corresponding He’s polynomials [ H1,n, H1,n ] [11] are H1,n = ∑∞ n=0 y1,iy1,n−i , H1,n = ∑∞ n=0 y1,iy1,n−i, n ≥ i, n = 0, 1, 2, . . . By solving the system (34) and (35), we obtain the iterations ỹi,0 (x, α), ỹi,1 (x, α) , . . . , ỹi,n (x, α) , i = 1, 2. Thus, the approximated solution in a series form is given by ỹi (x, α) = 2∑ n=0 ỹi,n (x, α) . W. Al-Hayani, M T. Younis / Eur. J. Pure Appl. Math, 16 (2) (2023), 1236-1259 1253 We compare the approximated solution using the HPM (n = 2) with the exact solution and the error produced in Tables 6–9. Table 10 demonstrates the maximum errors on the interval 0 ≤ x ≤ 1. Table 6: Comparison of the numerical results for y 1 (x, α) (Problem 2) α x y 1,E (x, α) y 1,2 (x, α) AEL 0.25 0.1 −0.0056156 −0.0056463 3.067E − 05 0.3 −0.0129290 −0.0128439 8.503E − 05 0.5 −0.0149820 −0.0146152 3.667E − 04 0.7 −0.0120790 −0.0113873 6.917E − 04 0.9 −0.0048957 −0.0042341 6.615E − 04 0.50 0.1 −0.0224625 −0.0226306 1.681E − 04 0.3 −0.0517160 −0.0519550 2.389E − 04 0.5 −0.0599281 −0.0598771 5.100E − 05 0.7 −0.0483163 −0.0477882 5.280E − 03 0.9 −0.0195831 −0.0189384 6.447E − 04 0.75 0.1 −0.0505406 −0.0509607 4.200E − 04 0.3 −0.1163610 −0.1171638 8.027E − 04 0.5 −0.1348384 −0.1353070 4.685E − 04 0.7 −0.1087117 −0.1084163 2.953E − 04 0.9 −0.0440621 −0.0433940 6.680E − 04 1.00 0.1 −0.0898500 −0.0906715 8.214E − 04 0.3 −0.2068641 −0.2085212 1.657E − 03 0.5 −0.2397127 −0.2409320 1.219E − 03 0.7 −0.1932653 −0.1932705 5.257E − 06 0.9 −0.0783326 −0.0775769 7.557E − 04 W. Al-Hayani, M T. Younis / Eur. J. Pure Appl. Math, 16 (2) (2023), 1236-1259 1254 Table 7: Comparison of the numerical results for y1 (x, α) (Problem 2) α x y1,E (x, α) y1,2 (x, α) AEU 0.25 0.1 −0.2246251 −0.2272738 2.648E − 03 0.3 −0.5171603 −0.5226018 5.441E − 03 0.5 −0.5992819 −0.6039739 4.691E − 03 0.7 −0.4831632 −0.4850771 1.913E − 03 0.9 −.01958317 −0.1952747 5.569E − 04 0.50 0.1 −0.1797001 −0.1816565 1.956E − 03 0.3 −0.4137282 −0.4176773 3.949E − 03 0.5 −0.4794255 −0.4825976 3.172E − 03 0.7 −0.3865306 −0.3873276 7.970E − 04 0.9 −0.1566653 −0.1556608 1.004E − 03 0.75 0.1 −0.1347751 −0.1361216 1.346E − 03 0.3 −0.3102962 −0.3129718 2.675E − 03 0.5 −0.3595691 −0.3615538 1.984E − 03 0.7 −0.2898979 −0.2900114 1.134E − 04 0.9 −0.1174990 −0.1163858 1.113E − 03 1.00 0.1 −0.0898500 −0.0906715 8.214E − 04 0.3 −0.2068641 −0.2085212 1.657E − 03 0.5 −0.2397127 −0.2409320 1.219E − 03 0.7 −0.1932653 −0.1932705 5.257E − 06 0.9 −0.0783326 −0.0775769 7.557E − 04 Table 8: Comparison of the numerical results for y 2 (x, α) (Problem 2) α x y 2,E (x, α) y 2,2 (x, α) AEL 0.25 0.1 0.0457031 0.0457084 5.327E − 06 0.3 0.1066406 0.1068146 1.740E − 04 0.5 0.1269531 0.1274189 4.657E − 04 0.7 0.1066406 0.1074269 7.862E − 04 0.9 0.0457031 0.0464363 7.332E − 04 0.50 0.1 0.0506250 0.0506059 1.900E − 05 0.3 0.1181250 0.1182502 1.252E − 04 0.5 0.1406250 0.1410804 4.554E − 04 0.7 0.1181250 0.1190405 9.155E − 04 0.9 0.0506250 0.0515624 9.374E − 04 0.75 0.1 0.0639843 0.0639516 3.277E − 05 0.3 0.1492968 0.1494366 1.397E − 04 0.5 0.1777343 0.1783339 5.995E − 04 0.7 0.1492968 0.1506143 1.317E − 03 0.9 0.0639843 0.0653802 1.395E − 03 1.00 0.1 0.0900000 0.0900165 1654E − 05 0.3 0.2100000 0.2103843 3843E − 04 0.5 0.2500000 0.2511554 1155E − 03 0.7 0.2100000 0.2122592 2259E − 03 0.9 0.0900000 0.0922510 2251E − 03 W. Al-Hayani, M T. Younis / Eur. J. Pure Appl. Math, 16 (2) (2023), 1236-1259 1255 Table 9: Comparison of the numerical results for y2 (x, α) (Problem 2) α x y2,E (x, α) y2,2 (x, α) AEU 0.25 0.1 0.1785937 0.1793249 7.311E − 04 0.3 0.4167187 0.4195254 2.806E − 03 0.5 0.4960937 0.5014287 5.335E − 03 0.7 0.4167187 0.4243124 7.593E − 03 0.9 0.1785937 0.1846205 6.026E − 03 0.50 0.1 0.16875 0.1692397 4.897E − 04 0.3 0.39375 0.3958039 2.053E − 03 0.5 0.46875 0.4728440 4.094E − 03 0.7 0.39375 0.3998311 6.081E − 03 0.9 0.16875 0.1737862 5.036E − 03 0.75 0.1 0.1420312 0.1422721 2.408E − 04 0.3 0.3314062 0.3326319 1.225E − 03 0.5 0.3945312 0.3972166 2.685E − 03 0.7 0.3314062 0.3357163 4.310E − 03 0.9 0.1420312 0.1458354 3.804E − 03 1.00 0.1 0.0900000 0.0900165 1654E − 05 0.3 0.2100000 0.2103843 3843E − 04 0.5 0.2500000 0.2511554 1155E − 03 0.7 0.2100000 0.2122592 2259E − 03 0.9 0.0900000 0.0922510 2251E − 03 Table 10: Maximum errors (Problem 2) ỹ1 (x, α) ỹ2 (x, α) α n y 1 (x, α) y1 (x, α) y 2 (x, α) y2 (x, α) 0.25 2 7.879E − 04 5.445E − 03 8.800E − 04 7.809E − 03 0.50 2 7.074E − 04 3.949E − 03 1.084E − 03 6.395E − 03 0.75 2 8.027E − 04 2.675E − 03 1.595E − 03 4.694E − 03 1.00 2 1.657E − 03 1.657E − 03 2.639E − 03 2.639E − 03 In Figures 9–12, a very good agreement is shown between the exact solution ( ỹi,E(x, 0.5), i = 1, 2 ) with a continuous line and the approximate solution by the HPM with Green’s function (ỹi,2(x, 0.5), i = 1, 2) with the symbol o. W. Al-Hayani, M T. Younis / Eur. J. Pure Appl. Math, 16 (2) (2023), 1236-1259 1256 Figure 9: Plot of y 1,E (x, 0.5) and y 1,2 (x, 0.5). Figure 10: Plot of Plot of y1,E(x, 0.5) and y1,2(x, 0.5). Figure 11: Plot of y 2,E (x, 0.5) and y 2,2 (x, 0.5). REFERENCES 1257 Figure 12: Plot of y2,E(x, 0.5) and y2,2(x, 0.5). 6. Conclusions The HPM with Green’s function has been used to solve linear and nonlinear fuzzy system of BVPs in this paper. The HPM we proved to be highly effective in address- ing problems we describe an efficient method for solving fuzzy system of boundary value problem. We employed the HPM technique using Green functions to simplify the compu- tation of the problem into linear and nonlinear fuzzy system of boundary value problem. This method is computationally efficient, as demonstrated by illustrative in two prob- lems. Moreover, we tested our approach through numerical simulations, and the results show that our method effectively solve the problems at hand. Overall, our findings sug- gest that the HPM approach with Green functions can be a valuable tool for addressing similar problems in the future. The results obtained when (α = 1) are compatible with those methods in the literature such as the Variational iteration method, Sinc-collocation method, and Boundary value methods. References [1] Ghazi Abed Meften, Ali Hasan Ali, Khalil S Al-Ghafri, Jan Awrejcewicz, and Omar Bazighifan. Nonlinear stability and linear instability of double-diffusive convection in a rotating with ltne effects and symmetric properties: brinkmann-forchheimer model. Symmetry, 14(3):565, 2022. [2] Mohammed J Ahmed and Waleed Al-Hayani. The homotopy perturbation method to solve initial value problems of first order with discontinuities. AL-Rafidain Journal of Computer Sciences and Mathematics, 16(2):61–70, 2022. [3] Rasha F. Ahmed, Waleed Al-Hayani, and Abbas Y. Al-Bayati. The homotopy analy- sis method to solve the nonlinear system of volterra integral equations and applying REFERENCES 1258 the genetic algorithm to enhance the solutions. European Journal of Pure and Ap- plied Mathematics, In Press, 2023. https://doi.org/10.29020/nybg.ejpam.v16i2. 4764. [4] M Tahmina Akter and MA Mansur Chowdhury. Homotopy perturbation method for solving nonlinear partial differential equations. IOSR Journal of Mathematics, 12(5):59–69, 2016. [5] Ali Hasan Ali, Ghazi Abed Meften, Omar Bazighifan, Mehak Iqbal, Sergio Elaskar, and Jan Awrejcewicz. A study of continuous dependence and symmetric properties of double diffusive convection: forchheimer model. Symmetry, 14(682):1–18, 2022. [6] Sarmad A Altaie, Nidal Anakira, Ali Jameel, Osama Ababneh, Ahmad Qazza, and Abdel Kareem Alomari. Homotopy analysis method analytical scheme for develop- ing a solution to partial differential equations in fuzzy environment. Fractal and Fractional, 6(8):419, 2022. [7] TA Biala and SN Jator. A family of boundary value methods for systems of second- order boundary value problems. International Journal of Differential Equations, 2017:1–12, 2017. [8] Snehashish Chakraverty, Nisha Mahato, Perumandla Karunakar, and Tha- rasi Dilleswar Rao. Advanced numerical and semi-analytical methods for differential equations. John Wiley & Sons, 2019. [9] Mohamed El-Gamel. Sinc-collocation method for solving linear and nonlinear system of second-order boundary value problems. Applied Mathematics, 3(11):1627–1633, 2012. [10] R Ezzati and S Ziari. Numerical solution and error estimation of fuzzy fredholm integral equation using fuzzy bernstein polynomials. Australian Journal of Basic and Applied Sciences, 5(9):2072–2082, 2011. [11] Asghar Ghorbani. Beyond adomian polynomials: he polynomials. Chaos, Solitons & Fractals, 39(3):1486–1492, 2009. [12] Ji-Huan He. A coupling method of a homotopy technique and a perturbation technique for non-linear problems. International journal of non-linear mechanics, 35(1):37–43, 2000. [13] Ji-Huan He. Homotopy perturbation method for solving boundary value problems. Physics letters A, 350(1-2):87–88, 2006. [14] Ji-Huan He and Yusry O El-Dib. Homotopy perturbation method with three expan- sions. Journal of Mathematical Chemistry, 59:1139–1150, 2021. [15] Shijun Liao. Beyond perturbation: introduction to the homotopy analysis method. CRC press, 2003. REFERENCES 1259 [16] Shijun Liao. Homotopy analysis method in nonlinear differential equations. Springer, 2012. [17] Junfeng Lu. Variational iteration method for solving a nonlinear system of second- order boundary value problems. Computers & Mathematics with Applications, 54(7- 8):1133–1138, 2007. [18] Ghazi Abed Meften and Ali Hasan Ali. Continuous dependence for double diffusive convection in a brinkman model with variable viscosity. Acta Universitatis Sapientiae, Mathematica, 14(1):125–146, 2022. [19] R Shanthi, T Iswarya, J Visuvasam, L Rajendran, and Michael EG Lyons. Voltam- metric and mathematical analysis of adsorption of enzymes at rotating disk electrode. Int. J. Electrochem. Sci, 17(220433):2, 2022. [20] Srinivasarao Thota. Solution of generalized abel’s integral equations by homotopy per- turbation method with adaptation in laplace transformation. Sohag J Math, 9(2):29– 35, 2022. [21] E Arul Vijayalakshmi, SS Santra, T Botmart, H Alotaibi, GB Loganathan, M Kan- nan, J Visuvasam, and V Govindan. Analysis of the magnetohydrodynamic flow in a porous medium. AIMS Mathematics, 7(8):15182–15194, 2022. [22] E Arul Vijayalakshmi, James Visuvasam, and M Kannan. Theoretical analysis of mhd mixed convection from a vertical plate embedded in a porous medium with a convective boundary condition. ECS Transactions, 107(1):18507–18522, 2022. [23] J Visuvasam, Angela Molina, E Laborda, and L Rajendran. Mathematical models of the infinite porous rotating disk electrode. Int. J. Electrochem. Sci, 13:9999–10022, 2018. [24] Mahasin Thabet Younis and Waleed Mohammed Al-Hayani. Solving fuzzy system of volterra integro-differential equations by using adomian decomposition method. European Journal of Pure and Applied Mathematics, 15(1):290–313, 2022.