EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 3, Article Number 6633 ISSN 1307-5543 – ejpam.com Published by New York Business Global A Boundary-Value Problem with Caputo-Hadamard Fractional Derivative: Analysis and Numerical Solution Afrah Sadiq Hasan1,∗, Shayma Adil Murad1 1 Department of Mathematics, College of Science, University of Duhok, Duhok, Iraq Abstract. We investigate a boundary-value problem governed by a fractional differential equa- tion, which is non-linear. The fractional derivative is the combined Caputo-Hadamard fractional derivative. We establish the required conditions for the existence and uniqueness of solutions, utilising the two standard fixed-point theorems, Banach fixed-point and Sadovskii fixed-point. Furthermore, we address and analyse the problem’s stability through demanding conditions using the Ulam-Hyers and Ulam-Hyers-Rassias stability methods. To illustrate the theoretical results, we present an example that validates the existence, uniqueness, and stability criteria. The analytical solution obtained from the problem is discretised with fractional rectangular, Lln,1 interpolation on a non-uniform mesh. The resulting system of non-linear equations is solved using the Newton- Raphson method, incorporating the Jacobian matrix to couple η(t) and θ(t, η(t)) non-linearly. The stability and reliability of the suggested numerical approach are examined through illustrative examples. 2020 Mathematics Subject Classifications: 34A08, 34K37, 34A12, 34Dxx, 65Lxx Key Words and Phrases: Lane-Emden fractional differential equation, existance, uniqueness, stability, numerical solution 1. Introduction It is an indisputable fact that fractional differential equations (FDEs) are crucial for modeling complex systems where traditional integer-order equations fall short, see [1, 2]. They offer a more flexible and accurate framework for governing processes with non- local behaviors and irregular dynamics. This makes FDEs valuable across various fields, including physics, biology, engineering, and finance, where systems exhibit non-linear and complex characteristics. The ability of FDEs to capture a wider range of phenomena makes them an essential tool for advancing our understanding of intricate real-world problems. See [3–8] and references therein. ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i3.6633 Email addresses: afrah.hasan@uod.ac (A. S. Hasan), Shayma.murad@uod.ac (S. A. Murad) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) A. S. Hasan, S. A. Murad / Eur. J. Pure Appl. Math, 18 (3) (2025), 6633 2 of 24 The classical Lane-Emden equation, rooted in astrophysics, was pioneered by [9] and [10] to model the equilibrium of self-gravitating polytropic gas spheres. Its dimensionless form, η′′(t) + µ t η′(t) = −ηn, 0 < t ≤ 1, µ > 0, (1) governs density profiles in stars, where η(t) represents scaled density and n is the polytropic index [11]. Chandrasekhar’s seminal work formalized its role in stellar structure theory, linking solutions (polytropes) to configurations of stars and gaseous planets. Prialnik in his book [12] further contextualized its applications in stellar evolution, particularly in modeling pre-main-sequence stars and degenerate cores. Beyond astrophysics, the equa- tion’s singular term t−1 and nonlinearity ηn inspired adaptations in plasma physics and radiative cooling. The Lane-Emden fractional differential equation (LEFDE) emerged as a modern extension, replacing integer derivatives with fractional operators to incorporate memory effects and anomalous transport. A generalized form, Dβ+αη(t) + µ t Dαη(t) = θ(t, η(t)), 0 < α, β < 1, (2) addresses non-local dynamics in systems like viscoelastic collapsing clouds or turbulent plasmas. Dβ+α and Dα are the fractional derivatives of order β + α and α, respectively. While [11], [13] and [12] focus on classical theory, recent studies leverage fractional calculus to resolve discrepancies in observational data, such as non-isothermal collapse in molec- ular clouds. This fractional framework retains the singular coefficient t−1 but introduces flexibility in modeling multi-scale phenomena, bridging gaps between classical polytropic assumptions and complex astrophysical systems. Lately, the study of initial and boundary-value problems that are governed by LEFDE has garnered substantial attention from many researchers. Ibrahim in [14] studied the existence of the non-linear LEFDE Dβ ( Dα + µ t ) η(t) = θ(t, η(t)), (3) for 0 < α, β ≤ 1, 0 < t ≤ 1, subject to the boundary conditions η(0) = η(t1) = η(1) = 0 for some t1 ∈ (0, 1). The same author in [15] studied the stability of the linear LEFDE Dβ ( Dα + µ t ) η(t) = θ(t), (4) for 0 < α, β ≤ 1, 0 < t ≤ 1, with boundary conditions η(0) = l1 and η(1) = l2, where l1 and l2 are constants. In [16] the authors used a combination of Chebyshev wavelets and a finite difference approaches to numerically solve the LEFDE Dαη(t) + µ tα−β Dβη(t) = θ(t, η(t)), (5) for 1 < α ≤ 2, 0 < β ≤ 1, 0 < t ≤ 1, subject to the initial or boundary conditions. In a study by [17], they considered the LEFDE in an n-dimensional system where each equation A. S. Hasan, S. A. Murad / Eur. J. Pure Appl. Math, 18 (3) (2025), 6633 3 of 24 consists of two arbitrary differential orders in terms of Caputo fractional derivative. They proved the existence and uniqueness using Krasnoselskii and Banach’s fixed point theo- rems and they showed that the system is both Ulam–Hyers and Ulam-Hyers-Rassias stable according to the proposed conditions. The existence and stability of LEFDE has been studied by [18] using advanced monotonicity, concentration-compactness, and Sobolev- type inequalities techniques. In [19], the author presents a modern analytical method for solving nonlinear singular type of LEFDE with Liouville–Caputo derivatives, focusing on the conditions that are proving the existence and uniqueness of solutions. The author combines techniques in fractional calculus alongside with advanced analytical methods to establish the conditions under which these solutions exist and are unique. The stability of solutions for two classes of LEFDE has been studied by [20] thought Lyapunov’s direct method while the existence and uniqueness of solutions are demonstrated using Banach’s fixed-point theory. The solution of LEFDE analytically studied by [21] using the tech- niques of power series approaches. Most recently, [22] explored the LEFDE with Caputo derivatives. The author managed to study the existence and uniqueness of mild solutions undergoing the Bielecki-type norm. In addition, the author manifested that, the mild solution is Ulam-Hyres type stable. Motivated by the above studies, we examine the existence, uniqueness, and stability of the fractional Lane-Emden Boundary-Value Problem. We consider the nonlinear Lane- Emden equation with multiple fractional derivatives: CHD β ( CHD α + µ t ) η(t) = θ(t, η(t)), (6) supplemented with boundary conditions: η(a) = ξ1, (7) η(1) = ξ2, (8) where µ, ξ1, ξ2 ∈ R, t ∈ [a, 1] and 0 < a < 1. In equation (6), α and β are positive non- integer numbers less than one; we obtain the classical Lane-Emden equation for α = β = 1. The two functions η(t) and θ(t, η(t)) belong to the set of all continuous function on [a, 1], i.e. η(t), θ(t, η(t)) ∈ C ([a, 1], R). Here CHD β represents the Caputo-Hadamard fractional derivative of order β as stated in Definition 3. We study the existence, uniqueness, and stability of the problem (6-8) on the interval [a, 1] such that 1/a is finite. In addition, we assemble a numerical scheme, by applying the fractional rectangular, Lln,1 interpolation on a logarithmic grid to approximate the derived analytical solution (18). The resulting system of nonlinear equation solved by the Newton-Raphson method using the Jacobian matrix to handle the non-linear coupling between η(t) and θ(t, η(t)). Numerical approaches to solve Caputo-Hadamard fractional differential equations are complicated due to the logarithmic kernel, which is weakly singular. Recently, some re- searchers are investigating the numerical schemes that can be considered. In a study A. S. Hasan, S. A. Murad / Eur. J. Pure Appl. Math, 18 (3) (2025), 6633 4 of 24 by [23], they examined finite difference methods for Caputo-Hadamard derivative frac- tional differential equations. In their study, the analogous Volterra integral equations were approximated via fractional rectangle, Llog,1 interpolation when they used the modified predictor–corrector for Caputo-Hadamard fractional differential equations. [24] used the local discontinuous Galerkin method to solve boundary-value problems with the Caputo- Hadamard fractional derivative. Furthermore, [25] numerically solved Caputo–Hadamard fractional differential equations with the graded meshes. The organization of this paper is as follows: Section 2 presents the essential definitions, theorems, remarks, and lemmas that support our current investigation. Section 3 addresses the existence and uniqueness of the solution for the problem (6-8). Section 4 examines the Ulam-Hyers and Ulam-Hyers-Rassias stability results, whilst Section 5 presents the nu- merical scheme using Newton-Raphson method to solve the Caputo-Hadamard fractional problem. Examples along with relevant graphs are presented in Sections 4-5. 2. Preliminaries This section presents key definitions, lemmas, and remarks that form the foundational concepts for this study and will be referenced and served as the groundwork for the analysis and discussions in this study. Definition 1. [7] For a given function η(t), the left-sided fractional integral of order β > 0 of Hadamard type is defined as HI β a η(t) = 1 Γ(β) ∫ t a ( ln ( t τ ))β−1 η(τ) τ dτ, (9) where Γ(·) refers to the Euler Gamma function. Definition 2. [7] For a given function η(t), the left-sided fractional derivative of order β > 0 (n− 1 < β < n ∈ N) of Hadamard type is defined as HD β aη(t) = 1 Γ(n− β) δn ∫ t a ( ln ( t τ ))n−β−1 η(τ) τ dτ, (10) where δn = ( t ddt )n . With changing the order of integration and derivative in (10) we arrive at the following. Definition 3. [26] If η(t) belongs to C([a, 1],R), then, the left-sided fractional derivative of order β > 0 of Caputo-Hadamard type is defined as CHD β aη(t) = 1 Γ(n− β) ∫ t a ( ln ( t τ ))n−β−1 δnη(τ) dτ τ , (11) where δn = ( τ d dτ )n and n− 1 < β < n. A. S. Hasan, S. A. Murad / Eur. J. Pure Appl. Math, 18 (3) (2025), 6633 5 of 24 Lemma 1. [27] For a given function η(t), if β > 0 (n− 1 < β < n) where n is an integer greater than 0, then, the Hadamard integral operates on the Caputo-Hadamard derivatives as HI β a+ CHD β a+ η(t) = η(t) + c0 + c1 ln ( t a ) + c2 ( ln ( t a ))2 + · · ·+ cn−1 ( ln ( t a ))n−1 , (12) where ci = δiη(a) i! , i = 0, 1, . . . , n− 1, are constants. Lemma 2. [27] For β > 0 and α > 0, we have HI β a+ ( ln ( t a ))α = Γ(α+ 1) Γ(α+ β + 1) ( ln ( t a ))α+β . (13) Theorem 1. [3] Let ϕ1 and ϕ2 be two opertors such that their combined mapping ϕ1+ϕ2 forms a k∗set-contraction for 0 ≤ k∗ < 1. This map is also condensing if the following conditions are satisfied: I- ϕ1, ϕ2 : B ⊆ Q → Q are operators on the Banach space Q where Q = C([a, 1],R). II- ϕ1 is k∗-contractive, that is, ∥ϕ1(x) − ϕ1(y)∥ ≤ k∗∥x − y∥ for all x and y in the domain and fixed k∗ ∈ [0, 1). III- ϕ2 is compact. Theorem 2. [3] (Sadovskii Fixed Point Theorem) For a subset B of a Banach space Q which is convex, bounded, and closed the condensing map ϕ : B → B has a fixed point. Theorem 3. [28] (Banach Fixed Point Theorem). For a continuous operator ϕ : B → B which is k∗-contractive there is a unique fixed point. Definition 4. [29] The boundary-value problem (6-8) is Ulam-Hyers stable if there exists a real constant ch > 0 such that for any ϵ > 0, and for every solution η(t) ∈ Q of the inequality ∣∣∣CHDβ ( CHD α + µ t ) η(t)− θ(t, η(t)) ∣∣∣ ≤ ϵ, (14) there exists a solution Z(t) ∈ Q of the problem (6-8) with |η(t)− Z(t)| ≤ chϵ, (15) for all t ∈ [a, 1]. Remark 1. [30] A function η(t) ∈ Q is a solution of the inequality (14) if and only if there exists a function g(t) ∈ Q (which is dependent on η), such that: (i) |g(t)| ≤ ϵ for all t ∈ [a, 1], A. S. Hasan, S. A. Murad / Eur. J. Pure Appl. Math, 18 (3) (2025), 6633 6 of 24 (ii) CHD β ( CHD α + µ t ) η(t) = θ(t, η(t)) + g(t). Definition 5. [29] The boundary-value problem (6-8) has the stability of the type Ulam- Hyers-Rassiass if there exists a real constant cλ > 0 such that for any ϵ > 0, and for every solution η(t) ∈ Q of the inequality∣∣∣CHDβ ( CHD α + µ t ) η(t)− θ(t, η(t)) ∣∣∣ ≤ ϵψ(t), (16) where ψ(t) ∈ Q is the control function, positive and non-decreasing, there exists a solution Z(t) ∈ Q of the problem (6-8) with |η(t)− Z(t)| ≤ cλϵψ(t), (17) for all t ∈ [a, 1]. Remark 2. [30] A function η(t) ∈ Q is a solution of the inequality (16) if and only if there exists a function g(t) ∈ Q (which is dependent on η), such that: (i) |g(t)| ≤ ϵψ(t) for all t ∈ [a, 1], (ii) CHD β ( CHD α + µ t ) η(t) = θ(t, η(t)) + g(t). Lemma 3. The solution of the boundary-value problem (6-8) is a function η(t) that belongs to C([a, 1],R) which takes the following form: η(t) = 1 Γ(β + α) ∫ t a ( ln ( t τ ))β+α−1 θ(τ, η(τ)) τ dτ − µ Γ(α) ∫ t a ( ln ( t τ ))α−1 η(τ) τ2 dτ − ( ln ( t a ))α bΓ(β + α) ∫ 1 a ( ln ( 1 τ ))β+α−1 θ(τ, η(τ)) τ dτ + µ ( ln ( t a ))α bΓ(α) ∫ 1 a ( ln ( 1 τ ))α−1 η(τ) τ2 dτ +Φ(t), (18) where Φ(t) = ξ1 + 1 b (ξ2 − ξ1) ( ln ( t a ))α and b = ( ln ( 1 a ))α . Proof By implementing the fractional Hadamard integral operator of order β, H a I β, given by Definition 1 on both sides of (6) and following Lemma 1 and Definition 1 we arrive at Dαη(t) = 1 Γ(β) ∫ t a ( ln ( t τ ))β−1 θ(τ, η(τ)) τ dτ − µ t η(t) + c̄, (19) where c̄ is a constant to be determined. Next, we operate both sides of (19) by the fractional Hadamard integral operator of order α, H a I α, through Lemma 1, to provide the following general solution A. S. Hasan, S. A. Murad / Eur. J. Pure Appl. Math, 18 (3) (2025), 6633 7 of 24 η(t) + c = 1 Γ(β + α) ∫ t a ( ln ( t τ ))β+α−1 θ(τ, η(τ)) τ dτ − µ Γ(α) ∫ t a ( ln ( t τ ))α−1 η(τ) τ2 dτ + c̄ Γ(α+ 1) ( ln ( t a ))α . (20) With applying the two boundary conditions (7-8) on equation (20) we obtain c = −ξ1 and c̄ = ξ2Γ(α+ 1) b − ξ1Γ(α+ 1) b − Γ(α+ 1) bΓ(β + α) ∫ 1 a ( ln ( 1 τ ))β+α−1 θ(τ, η(τ)) τ dτ + αµ b ∫ 1 a ( ln ( 1 τ ))α−1 η(τ) τ2 dτ. (21) After plugging in the known values of c̄ and c into (20), it turns out that the solution to the boundary-value problem (6-8) is given by (18). Through direct computation, the converse is obtained. The proof for Lemma 3 has been finalised. 3. Existence and Uniqueness of Solution In this section, first we apply the Banach fixed point theorem to show the existence and uniqueness of the boundary-value problem (6-8), and then, the Sadovoskii fixed point theorem for the existence of the solution to our problem. Let Br = {η(t) ∈ Q : ∥η∥ ≤ r} be a closed, bounded, and convex subset of Q = C([a, 1],R) where r is a positive constant such that r ≥ 2bΩ Γ(β+α+1) ( ln ( 1 a ))β + ∥Φ∥ 1− 2bµ aΓ(α+1) , 2bµ < aΓ(α+ 1). (22) Here Q is a Banach space of all continuous function from [a, 1] to R with the norm ∥η∥ = sup t {|η(t)|, t ∈ [a, 1]} for all η ∈ Q. (23) We proceed with Section 3 under the following two assumptions, namely H1 and H2, which state: H1: There exists a constant Ω > 0, such that |θ(t, η(t))| ≤ Ω for all t ∈ [a, 1] and η ∈ Q. H2: There exists a constant k > 0, such that ∥θ(t, η1) − θ(t, η2)∥ ≤ k∥η1 − η2∥ for all t ∈ [a, 1] and η1, η2 ∈ Q if: 0 < k Γ(β + α+ 1) ( ln ( 1 a ))β + µ aΓ(α+ 1) < 1 2b . We initiate by calling the operator ϕ : Q −→ Q to be set as: ϕ (η(t)) = 1 Γ(β + α) ∫ t a ( ln ( t τ ))β+α−1 θ(τ, η(τ)) τ dτ A. S. Hasan, S. A. Murad / Eur. J. Pure Appl. Math, 18 (3) (2025), 6633 8 of 24 − µ Γ(α) ∫ t a ( ln ( t τ ))α−1 η(τ) τ2 dτ − ( ln ( t a ))α bΓ(β + α) ∫ 1 a ( ln ( 1 τ ))β+α−1 θ(τ, η(τ)) τ dτ + µ ( ln ( t a ))α bΓ(α) ∫ 1 a ( ln ( 1 τ ))α−1 η(τ) τ2 dτ +Φ(t), t ∈ [a, 1]. (24) It has to be shown that ϕ has a fixed point. The fixed point is a solution to the boundary- value problem (6-8). We have to manifest that ϕ (Br) ⊂ Br. For any η ∈ Br: ∥ϕ (η) ∥ = Sup t ∣∣∣∣ 1 Γ(β + α) ∫ t a ( ln ( t τ ))β+α−1 θ(τ, η(τ)) τ dτ − µ Γ(α) ∫ t a ( ln ( t τ ))α−1 η(τ) τ2 dτ − ( ln ( t a ))α bΓ(β + α) ∫ 1 a ( ln ( 1 τ ))β+α−1 θ(τ, η(τ)) τ dτ + µ ( ln ( t a ))α bΓ(α) ∫ 1 a ( ln ( 1 τ ))α−1 η(τ) τ2 dτ +Φ(t) ∣∣∣∣, t ∈ [a, 1]. (25) Utilizing the assumption H1 and the condition (22), we arrive at: ∥ϕ (η) ∥ ≤ Ωφ1 + rφ2 + ∥Φ∥ ≤ r, (26) where φ1 = 2b Γ(β+α+1) ( ln ( 1 a ))β and φ2 = 2bµ aΓ(α+1) , and inequality (26) reveal that ϕ (Br) ⊂ Br. Progressing from the prior, we now prove the contraction requirement on the mapping. For any two functions η1 and η2 in Br, the norm of their difference follows: ∥ϕ (η2)− ϕ (η1)∥ = ∥∥∥∥∥ 1 Γ(β + α) ∫ t a ( ln ( t τ ))β+α−1 [θ(τ, η2(τ))− θ(τ, η1(τ))] dτ τ − µ Γ(α) ∫ t a ( ln ( t τ ))α−1 [η2(τ)− η1(τ)] dτ τ2 − ( ln ( t a ))α bΓ(β + α) ∫ 1 a ( ln ( 1 τ ))β+α−1 [θ(τ, η2(τ))− θ(τ, η1(τ))] dτ τ + µ ( ln ( t a ))α bΓ(α) ∫ 1 a ( ln ( 1 τ ))α−1 [η2(τ)− η1(τ)] dτ τ2 ∥∥∥∥∥ , t ∈ [a, 1]. (27) Equation (27) follows: ∥ϕ (η2)− ϕ (η1)∥ ≤ ( ln ( t a ))α+β Γ(β + α+ 1) ∥θ2 − θ1∥+ µ ( ln ( t a ))α aΓ(α+ 1) ∥η2 − η1∥ A. S. Hasan, S. A. Murad / Eur. J. Pure Appl. Math, 18 (3) (2025), 6633 9 of 24 + ( ln ( t a ))α ( ln ( 1 a ))α+β bΓ(β + α+ 1) ∥θ2 − θ1∥+ µ ( ln ( t a ))α ( ln ( 1 a ))α abΓ(α+ 1) ∥η2 − η1∥ . = φ1 ∥θ2 − θ1∥+ φ2 ∥η2 − η1∥ , (28) Using assumption (H2), (28) arrives at: ∥ϕ (η2)− ϕ (η1)∥ ≤ Υ ∥η2 − η1∥ , (29) where 0 < Υ = kφ1 + φ2 < 1. Using Theorem 3, we can show that the solution to the boundary-value problem (6–8) is unique. For the existence of a solution for the boundary- value problem (6-8), we apply the Sadovoskii fixed point theorem. We consider the oper- ator to be ϕ, as it is defined in 24. At this point, we introduce two operators called ϕ1 and ϕ2 that map Q to itself as follows: ϕ1 (η(t)) =− µ Γ(α) ∫ t a ( ln ( t τ ))α−1 η(τ) τ2 dτ + µ ( ln ( t a ))α bΓ(α) ∫ 1 a ( ln ( 1 τ ))α−1 η(τ) τ2 dτ +Φ(t), t ∈ [a, 1], (30) and ϕ2 (η(t)) = 1 Γ(β + α) ∫ t a ( ln ( t τ ))β+α−1 θ(τ, η(τ)) τ dτ − ( ln ( t a ))α bΓ(β + α) ∫ 1 a ( ln ( 1 τ ))β+α−1 θ(τ, η(τ)) τ dτ, t ∈ [a, 1], (31) respectively. Instead of seeking the existence of a fixed point for the operator ϕ (η(t)), we shall do the same through the sum of the two previously defined operators in (30-31). We will use Theorem 2 and the contraction condition 0 ≤ µ(b + 1)/(aΓ(α + 1)) < 1 to show that ϕ1 + ϕ2 has a fixed point using the steps below. Step 1: It has previously manifested that ϕ (Br) ⊂ Br. Step 2: We are compelled to exhibit that ϕ2 is compact. For any t1, t2 ∈ [a, 1] with t1 < t2 and for any η(t) ∈ Br we have: |ϕ2(η(t2))− ϕ2(η(t1))| = ∣∣∣∣∣ 1 Γ(β + α) ∫ t2 a ( ln ( t2 τ ))β+α−1 θ(τ, η(τ)) τ dτ − ( ln ( t2 a ))α bΓ(β + α) ∫ 1 a ( ln ( 1 τ ))β+α−1 θ(τ, η(τ)) τ dτ − 1 Γ(β + α) ∫ t1 a ( ln ( t1 τ ))β+α−1 θ(τ, η(τ)) τ dτ + ( ln ( t1 a ))α bΓ(β + α) ∫ 1 a ( ln ( 1 τ ))β+α−1 θ(τ, η(τ)) τ dτ ∣∣∣∣∣ . (32) A. S. Hasan, S. A. Murad / Eur. J. Pure Appl. Math, 18 (3) (2025), 6633 10 of 24 ≤ 1 Γ(β + α) ∫ t1 a [( ln ( t2 τ ))β+α−1 − ( ln ( t1 τ ))β+α−1 ] |θ(τ, η(τ))| τ dτ + 1 Γ(β + α) ∫ t2 t1 ( ln ( t2 τ ))β+α−1 |θ(τ, η(τ))| τ dτ + 1 bΓ(β + α) ∣∣∣∣(ln( t1a ))α − ( ln ( t2 a ))α∣∣∣∣ ∫ 1 a ( ln ( 1 τ ))β+α−1 |θ(τ, η(τ))| τ dτ. (33) It is shown in (33) that as |t2 − t1| −→ 0 we have |ϕ2(η(t2))− ϕ1(η(t1))| −→ 0. Hence ϕ is equicontinuous. Since ϕ (Br) ⊂ Br, then, ϕ2 is uniformly bounded. As a consequence of the Arzelà-Ascoli theorem, ϕ2 (Br) is compact. Step 3: For this part, we have to demonstrate that ϕ1 is k∗-contractive. For any η1 and η2 in Br, we have: ∥ϕ1 (η1)− ϕ1 (η2)∥ ≤ Sup t { µ Γ(α) ∫ t a ( ln ( t τ ))α−1 |η1(τ)− η2(τ)| τ2 dτ + µ ( ln ( t a ))α bΓ(α) ∫ 1 a ( ln ( 1 τ ))α−1 |η1(τ)− η2(τ)| τ2 dτ } , t ∈ [a, 1]. ≤ µb aΓ(α+ 1) ∥η1 − η2∥+ µ aΓ(α+ 1) ∥η1 − η2∥ = k∗ ∥η1 − η2∥ . (34) Where k∗ = µ(b+1) aΓ(α+1) . Step 4: In the last part, we need to show that ϕ is condensing. Since ϕ1 is continuous and k∗-contractive, and we found that ϕ2 is compact, therefore, by Theorem 1, we observe that ϕ = ϕ1 + ϕ2 is a condensing map on Br. From the above four steps and by the Sadovski Theorem 2, we arrive at the conclusion that the map ϕ has a fixed point. 4. Stability This section examines the stability of the boundary-value problem (6-8) utilising two established definitions of stability: Ulam-Heyrs stability and Ulam-Heyrs-Rassias stability. Theorem 4. Assume that θ : [a, 1]×R −→ R is a continuous function and the assumption H2 holds. Then, the solution of the boundary-value problem (6-8) is Ulam-Hyers stable. Proof: Let η(t) ∈ C([a, 1],R) be a solution of the inequality (14) which satisfies boundary conditions (7-8). Through Remark 1, we have CHD β ( Dα + µ t ) η(t) = θ(t, η(t)) + g(t), (35) A. S. Hasan, S. A. Murad / Eur. J. Pure Appl. Math, 18 (3) (2025), 6633 11 of 24 where g(t) possesses the same property mentioned in Remark 1. The solution of the perturbed problem (35) supplemented with (7-8) can be found to be: η(t) = 1 Γ(β + α) ∫ t a ( ln ( t τ ))β+α−1 θ(τ, η(τ)) τ dτ − µ Γ(α) ∫ t a ( ln ( t τ ))α−1 η(τ) τ2 dτ + 1 Γ(β + α) ∫ t a ( ln ( t τ ))β+α−1 g(τ) τ dτ − ( ln ( t a ))α bΓ(β + α) ∫ 1 a ( ln ( 1 τ ))β+α−1 θ(τ, η(τ)) τ dτ + µ ( ln ( t a ))α bΓ(α) ∫ 1 a ( ln ( 1 τ ))α−1 η(τ) τ2 dτ − ( ln ( t a ))α bΓ(β + α) ∫ 1 a ( ln ( 1 τ ))β+α−1 g(τ) τ dτ +Φ(t), (36) where (36) satisfies the following inequality:∣∣∣∣∣η(t)− 1 Γ(β + α) ∫ t a ( ln ( t τ ))β+α−1 θ(τ, η(τ)) τ dτ + µ Γ(α) ∫ t a ( ln ( t τ ))α−1 η(τ) τ2 dτ + ( ln ( t a ))α bΓ(β + α) ∫ 1 a ( ln ( 1 τ ))β+α−1 θ(τ, η(τ)) τ dτ − µ ( ln ( t a ))α bΓ(α) ∫ 1 a ( ln ( 1 τ ))α−1 η(τ) τ2 dτ − Φ(t) ∣∣∣∣∣ ≤ 2 ϵ b ( ln ( 1 a ))β Γ(β + α+ 1) , (37) for all t ∈ [a, 1]. Let Z(t) be the unique solution for our boundary-value problem (6-8). The left side of inequality (15) conforms to |η(t)− Z(t)| ≤ 2 ϵ ( ln ( 1 a ))β+α Γ(β + α+ 1) + 1 Γ(β + α) ∫ t a ( ln ( t τ ))β+α−1 |θ(τ, η(τ))− θ(τ, Z(τ))| dτ τ + µ Γ(α) ∫ t a ( ln ( t τ ))α−1 |η(τ)− Z(τ)| dτ τ2 + ( ln ( t a ))α bΓ(β + α) ∫ 1 a ( ln ( 1 τ ))β+α−1 |θ(τ, η(τ))− θ(τ, Z(τ))| dτ τ + µ ( ln ( t a ))α bΓ(α) ∫ 1 a ( ln ( 1 τ ))α−1 |η(τ)− Z(τ)| dτ τ2 . (38) Considering the Lipchitz condition H2, the inequality (38) reads A. S. Hasan, S. A. Murad / Eur. J. Pure Appl. Math, 18 (3) (2025), 6633 12 of 24 |η(t)− Z(t)| ≤ 2bϵ ( ln ( 1 a ))β Γ(β + α+ 1) + (kφ1 + φ2) |η(t)− Z(t)| , (39) and by rewriting inequality (39), we obtain |η(t)− Z(t)| ≤ chϵ, (40) where ch = φ1 1−Υ . Inequality (40) via Definition 4 ensures that the boundary-value problem (6-8) is stable of type Ulam-Hayrs. To investigate the Ulam-Hayrs-Rassias type stability, we imposed a new assumption on the problem (6-8). H3: The control function ψ(t) in Definition 5 is a positive non-decreasing function such that: H a I α+βψ(t) ≤ λψ(t), (41) where λ > 0. Theorem 5. Assume that θ : [a, 1]×R −→ R is a continuous function and both hypotheses H2 and H3 hold. Then, the solution of the boundary-value problem (6-8) is Ulam-Hyers- Rassias type stable. Proof: Let η(t) ∈ C([a, 1],R) be a solution of the inequality (16) which satisfies boundary conditions (7-8). Through Remark 2, we obtain the perturbed problem (35), where |g(t)| is bounded by ϵψ(t) as stated in Remark 2, and ψ(t) is the control function that fulfills assumption H3. Thus we have∣∣∣∣∣η(t)− 1 Γ(β + α) ∫ t a ( ln ( t τ ))β+α−1 θ(τ, η(τ)) τ dτ + µ Γ(α) ∫ t a ( ln ( t τ ))α−1 η(τ) τ2 dτ + ( ln ( t a ))α bΓ(β + α) ∫ 1 a ( ln ( 1 τ ))β+α−1 θ(τ, η(τ)) τ dτ − µ ( ln ( t a ))α bΓ(α) ∫ 1 a ( ln ( 1 τ ))α−1 η(τ) τ2 dτ − Φ(t) ∣∣∣∣∣ ≤ 2ϵλψ(t), (42) for all t ∈ [a, 1]. If we call Z(t) the unique solution for our boundary-value problem (6-8), then, the left side of inequality (17) conforms to |η(t)− Z(t)| ≤ 2ϵλψ(t) + 1 Γ(β + α) ∫ t a ( ln ( t τ ))β+α−1 |θ(τ, η(τ))− θ(τ, Z(τ))| dτ τ + µ Γ(α) ∫ t a ( ln ( t τ ))α−1 |η(τ)− Z(τ)| dτ τ2 A. S. Hasan, S. A. Murad / Eur. J. Pure Appl. Math, 18 (3) (2025), 6633 13 of 24 + ( ln ( t a ))α bΓ(β + α) ∫ 1 a ( ln ( 1 τ ))β+α−1 |θ(τ, η(τ))− θ(τ, Z(τ))| dτ τ + µ ( ln ( t a ))α bΓ(α) ∫ 1 a ( ln ( 1 τ ))α−1 |η(τ)− Z(τ)| dτ τ2 . (43) Considering the second assumption H2, the inequality (43) reads |η(t)− Z(t)| ≤ 2ϵλψ(t) + (kφ1 + φ2) |η(t)− Z(t)| , (44) and by rewriting inequality (44), we obtain |η(t)− Z(t)| ≤ ϵcλψ(t), (45) where cλ = 2λ (1−Υ) . Inequality (45) guarantees the stability of type Ulam-Hyers-Rassias according to Definition 5. Example 1: Consider the following multifractional boundary-value problem: CHD β ( CHD α + µ t ) η(t) = |η(t)| e−2t (10 + t2) (1 + η(t)2) . (46) η(a) = 0, (47) η(1) = 0.1. (48) In this problem, (46-48), which is defined on [a, 1] where 0 < a < 1, we have θ(t, η(t)) = |η(t)|e−2t (10+t2)(1+η2) which satisfies both hypotheses H1 and H2 with Ω = k = e−2a 10+a2 , i.e.: |θ(t, η(t))| ≤ Ω = e−2a 10 + a2 , (49) and ∥θ(t, η1)− θ(t, η2)∥ ≤ e−2a 10 + a2 ∥η1 − η2∥. (50) Hence by Banach’s Fixed Point Theorem there exists a unique solution to the problem (46-48). To analyze the problem (46-48), we set the interval [a, 1] with a = 0.8, and a = 0.9, which read k ≈ 0.01898, and k ≈ 0.01529, respectively.Also, we set µ = 0.1. Next, we take the values of α and β to run on [0.1, 0.9]. We demonstrate the proposed contraction condition in inequality (29), Υ = kφ1 + φ2, which depends on α, β, µ, and a. A. S. Hasan, S. A. Murad / Eur. J. Pure Appl. Math, 18 (3) (2025), 6633 14 of 24 (a) (b) Figure 1: The range of contraction parameter Υ, with µ = 0.1 as α and β move from 0.1 to 0.9 for the values of a: (a) a = 0.8, and (b) a = 0.9. For stability, we assume that η(t) and Z(t) are the solutions for the perturbed (for some ϵ > 0) and unperturbed problem (46-48), respectively. Hence, η(t) is a solution for the inequality (14). Since both hypotheses H1 and H2 are satisfied through (49-50), then, by Theorem 4, both η(t) and Z(t) satisfy inequality (15), and consequently the problem (46-48) is Ulam-Hyers stable with ch = φ1 1−Υ . A. S. Hasan, S. A. Murad / Eur. J. Pure Appl. Math, 18 (3) (2025), 6633 15 of 24 The last part is to establish how the problem is stable by the means of Ulam-Hyers- Rassiass as stated in Definition 5. If η(t) is a solution to the inequality (16) for some constant ϵ > 0 and the control function is ψ(t) = ln ( t a ) , for t ∈ [a, 1], then, ψ(t) satisfies the hypothesis H3. In other words, H a I α+βψ(t) = H a I α+β ln ( t a ) (51) = Γ(2) Γ(β + α+ 2) ( ln ( t a ))α+β+1 = 1 Γ(β + α+ 2) ( ln ( t a ))α+β ln ( t a ) ≤ λψ(t), where λ = γ Γ(β+α+2) and γ represents the maximum value of ( ln ( t a ))α+β on [a, 1]. The value of λ depends on three parameters, a, α, and β. For example, for a = 0.9, α = β = 0.5, we have λ = 0.05268. Consequently, by Theorem 5, the problem (46-48) is Ulam-Hyers- Rassiass stable as the following inequality holds, |η(t)− Z(t)| ≤ ϵcλψ(t), (52) where cλ = 2λ (1−kφ1−φ2) . Both parameters, ch and cλ are depicted in Figure 2 for both a = 0.9 and a = 0.8. A. S. Hasan, S. A. Murad / Eur. J. Pure Appl. Math, 18 (3) (2025), 6633 16 of 24 (a) (b) Figure 2: The range of values of ch and cλ as both orders, α and β, grow from 0.1 up to 0.9 for a = 0.9 (left) and a = 0.8 (right), (a) ch, and (b) cλ. 5. Numerical Scheme In this section we introduce a numerical scheme to approximately solve the problem (6-8) through the Newton-Raphson method. We apply the fractional rectangular, Lln,1, interpolation on a logarithmic grid to approximate the derived analytical solution (18). A. S. Hasan, S. A. Murad / Eur. J. Pure Appl. Math, 18 (3) (2025), 6633 17 of 24 This results in a system of non-linear equations which will be solved by the Newton- Raphson method using the Jacobian matrix to handle the non-linear coupling between η(t) and θ(t, η(t)). We subdivide the domain using a non-uniform mesh on [a, 1], B = {t0, t1, · · · , tN}, for some positive integer N where t0 = a and tN = 1. The following formula is used to generate the graded mesh: ln(ti) = ln(t0) + ∆ ( i N )s , (53) where ∆ = ln(tN ) − ln(t0), for the considered interval we have ∆ = −ln(t0), and s is the mesh graded parameter which is sensitive. For s = 1 we obtain the uniform mesh (in the logarithmic scale), while for optimized s > 1, both the resolution of the non-local effects inherent to fractional derivatives and the logarithmic kernel close to t = a and the convergence rate are improved. For i = 0, formula (53) generates ln(t0) and for i = N , it generates ln(tN ). The numerical solution of (18) is equivalent to solving the two-point fractional boundary- value problem (6-8). At t = tq, (0 ≤ q ≤ N), tq ∈ B, denote ηq ≈ η(tq), we have: ηq = 1 Γ(β + α) ∫ tq a ( ln ( tq τ ))β+α−1 θ(τ, η(τ)) τ dτ − µ Γ(α) ∫ tq a ( ln ( tq τ ))α−1 η(τ) τ2 dτ − ( ln ( tq a ))α bΓ(β + α) ∫ 1 a ( ln ( 1 τ ))β+α−1 θ(τ, η(τ)) τ dτ + µ ( ln ( tq a ))α bΓ(α) ∫ 1 a ( ln ( 1 τ ))α−1 η(τ) τ2 dτ +Φq, (54) where Φq = Φ(tq) = ξ1 + 1 b (ξ2 − ξ1) ( ln ( tq a ))α . Furthermore, ηq = 1 Γ(β + α) q∑ i=1 ∫ ti ti−1 ( ln ( tq τ ))β+α−1 θ(τ, η(τ)) dτ τ − µ Γ(α) q∑ i=1 ∫ ti ti−1 ( ln ( tq τ ))α−1 η(τ) τ dτ τ + f(tq), (55) where f(tq) = − ( ln ( tq a ))α bΓ(β + α) N∑ i=1 ∫ ti ti−1 ( ln ( 1 τ ))β+α−1 θ(τ, η(τ)) dτ τ + µ ( ln ( tq a ))α bΓ(α) N∑ i=1 ∫ ti ti−1 ( ln ( 1 τ ))α−1 η(τ) τ dτ τ +Φq. (56) A. S. Hasan, S. A. Murad / Eur. J. Pure Appl. Math, 18 (3) (2025), 6633 18 of 24 Note that, the first sum of (56) is the same as the first sum of (55), for q = N . The same is true for the second sums of (55) and (56). Applying the fractional rectangular, Lln,1 interpolation to approximate the nonlinear terms θ(τ, η(τ)) and η(τ) τ , that is θ(τ, η(τ)) ≈ ln ( τ ti ) ln ( ti−1 ti )θi−1 + ln ( τ ti−1 ) ln ( ti ti−1 )θi, (57) η(τ) ≈ ln ( τ ti ) ln ( ti−1 ti )ηi−1 + ln ( τ ti−1 ) ln ( ti ti−1 )ηi, (58) where θi = θ(ti, η(ti)) and ηi = η(ti). Now, in conjunction with (57-58) and ηq in (55) we arrive at: ηq = 1 Γ(β + α) q∑ i=1 ∫ ti ti−1 ( ln ( tq τ ))β+α−1  ln ( τ ti ) ln ( ti−1 ti )θi−1 + ln ( τ ti−1 ) ln ( ti ti−1 )θi  dτ τ − µ Γ(α) q∑ i=1 ∫ ti ti−1 ( ln ( tq τ ))α−1  ln ( τ ti ) ln ( ti−1 ti )ηi−1 + ln ( τ ti−1 ) ln ( ti ti−1 )ηi  dτ τ2 + f(tq) = 1 Γ(β + α) q∑ i=1  θi−1 ln ( ti−1 ti ) ∫ ti ti−1 ( ln ( tq τ ))β+α−1 ln ( τ ti ) dτ τ + θi ln ( ti ti−1 ) ∫ ti ti−1 ( ln ( tq τ ))β+α−1 ln ( τ ti−1 ) dτ τ  − µ Γ(α) q∑ i=1  ηi−1 ln ( ti−1 ti ) ∫ ti ti−1 ( ln ( tq τ ))α−1 ln ( τ ti ) dτ τ2 + ηi ln ( ti ti−1 ) ∫ ti ti−1 ( ln ( tq τ ))α−1 ln ( τ ti−1 ) dτ τ2 + f(tq), (59) proceeding the integration in (59) using the change of variables tq = τeu we obtain: ηq = 1 Γ(β + α) q∑ i=1  θi−1 ln ( ti−1 ti ) ln ( tq ti ) (ln tq ti−1 )β+α − ( ln tq ti )β+α β + α − θi−1 ln ( ti−1 ti ) (ln tq ti−1 )β+α+1 − ( ln tq ti )β+α+1 β + α+ 1 A. S. Hasan, S. A. Murad / Eur. J. Pure Appl. Math, 18 (3) (2025), 6633 19 of 24 + θi ln ( ti ti−1 ) ln ( tq ti−1 ) (ln tq ti−1 )β+α − ( ln tq ti )β+α β + α − θi ln ( ti ti−1 ) (ln tq ti−1 )β+α+1 − ( ln tq ti )β+α+1 β + α+ 1  − µ Γ(α) q∑ i=1  (−1)αηi−1 tq ln ( ti−1 ti )[Γ(α+ 1,− ln ( tq ti−1 )) − Γ ( α+ 1,− ln ( tq ti ))] −(−1)α−1ηi−1 tq ln ( ti−1 ti ) ln ( tq ti−1 )[ Γ ( α,− ln ( tq ti−1 )) − Γ ( α,− ln ( tq ti ))] + (−1)α−1ηi tq ln ( ti ti−1 ) ln ( tq ti )[ Γ ( α,− ln ( tq ti−1 )) − Γ ( α,− ln ( tq ti ))] − (−1)αηi tq ln ( ti ti−1 ) [Γ(α+ 1,− ln ( tq ti−1 )) − Γ ( α+ 1,− ln ( tq ti ))]] , (60) where Γ(., .) is the incomplete gamma function. Rearranging (60) ηq + µ Γ(α) q∑ i=1 [νqiηi−1 + ν̄qiηi] = 1 Γ(β + α) q∑ i=1 [ωqiθi−1 + ω̄qiθi] + f(tq), (61) where νqi = (−1)α tq ln ( ti−1 ti )[Γ(α+ 1,− ln ( tq ti−1 )) − Γ ( α+ 1,− ln ( tq ti ))] − (−1)α−1 tq ln ( ti−1 ti ) ln ( tq ti−1 )[ Γ ( α,− ln ( tq ti−1 )) − Γ ( α,− ln ( tq ti ))] , ν̄qi = (−1)α−1 tq ln ( ti ti−1 ) ln ( tq ti )[ Γ ( α,− ln ( tq ti−1 )) − Γ ( α,− ln ( tq ti ))] − (−1)α tq ln ( ti ti−1 )[Γ(α+ 1,− ln ( tq ti−1 )) − Γ ( α+ 1,− ln ( tq ti ))] , ωqi = ln ( tq ti ) (ln tq ti−1 )β+α − ( ln tq ti )β+α (β + α) ln ( ti−1 ti ) − ( ln tq ti−1 )β+α+1 − ( ln tq ti )β+α+1 (β + α+ 1) ln ( ti−1 ti ) , ω̄qi = ln ( tq ti−1 ) (ln tq ti−1 )β+α − ( ln tq ti )β+α (β + α) ln ( ti ti−1 ) − ( ln tq ti−1 )β+α+1 − ( ln tq ti )β+α+1 (β + α+ 1) ln ( ti ti−1 ) . (62) A. S. Hasan, S. A. Murad / Eur. J. Pure Appl. Math, 18 (3) (2025), 6633 20 of 24 Based on the right-hand side of (61), we introduce the vector H = [η0, η1, . . . , ηN ]T . We further introduce an (N + 1) × (N + 1) matrix A. The first and last rows of the matrix A are [1, 0, 0, . . . , 0] and [0, 0, . . . , 0, 1], respectively, as a result of enforcing the boundary conditions (7) and (8) on η0 and ηN , respectively. The elements of the interior rows are represented by (Aqi) N i=1 with q = 2, . . . , N − 1, so that: Aqi = µ Γ(α)  νqi ; i = 1, νqi + ν̄qi ; i = 2, 3, . . . , q − 1, Γ(α)µ−1 + ν̄qi ; i = q, 0 ; q < i ≤ N. (63) For the left-hand side of (61), we assemble the vector B = [B0, B1, . . . , BN ]T , where Bq =  ξ1 ; q = 0, 1 Γ(β+α) ∑q i=1 [ωqiθi−1 + ω̄qiθi] + f(tq) ; q = 1, 2, . . . , N − 1, ξ2 ; q = N. (64) Therefore, the two-point boundary-value problem (6-8) numerically will be solved by im- plementing the Newton-Raphson iteration on the non-linear system: AH = B. (65) We rearrange the nonlinear system in (65) as: F(H) = AH− B = 0, (66) and we use a Matlab code to solve iteratively the following: H(k+1) = H(k) − J(H(k)) −1 F(H(k)). (67) Here, J(H) is the Jacobian matrix with entries Jpi = Api− ∂Bq ∂ηi . Initially, we start with H(0), which is generated by applying the linear interpolation ηp = ξ1 + (tp − a)/(1− a)(ξ2 − ξ1) using the boundary conditions (7) and (8). Then, from (67) we obtain H(1) and so on. The iteration stops when we achieve the required tolerance. In the following examples we set s to 2−β β (see [31]). Example 2: Consider the following multifractional boundary-value problem: CHD β ( CHD α + µ t ) η(t) = η(t) t ( 2 + µ erf (√ ln ( η(t) at ))) , (68) η(a) = a2, (69) η(1) = 1. (70) A. S. Hasan, S. A. Murad / Eur. J. Pure Appl. Math, 18 (3) (2025), 6633 21 of 24 Here, the non-linear function θ(t, η(t)) = η(t) t ( 2 + µ erf (√ ln η(t) at )) , where erf(·) is the error function. For α = β = 0.5 and µ = 0.1, η(t) = t2 is the exact solution for the problem (68-70). The exact solution of the problem (68-70) versus the numerical solution is depicted in Figures (3)a-b. Both graphs presented with N = 512, and the convergence was achieved with tolerance 1e-06 at the 4th iteration. (a) (b) Figure 3: Exact solution, η(t) = t2, (solid line) versus numerical solution (dashed line) for the problem (68-70), with µ = 0.1, α = β = 0.5, N = 512, (a) a = 0.9, and (b) a = 0.8. The reduction of the absolute error is depicted in Figures (4)a-b, which is of order O(1e-03), as the mesh increases. (a) (b) Figure 4: Absolute error for the problem (68-70) as the mesh number (N) increases: (a) for a = 0.9 at t = 0.92, and (b) for a = 0.8 at t = 8.2. Example 3: Consider the following multifractional boundary-value problem: CHD β ( CHD α + µ t ) η(t) = η(t) t + 2µ√ πt √ η(t) + 1. (71) η(a) = 0, (72) η(1) = ln ( 1 a ) . (73) Here, the non-linear function θ(t, η(t)) = η(t) t + 2µ√ πt √ η(t) + 1. For α = β = 0.5 and µ = 0.1, the exact solution for the problem (71-73) is η(t) = t ln ( t a ) . The exact solution A. S. Hasan, S. A. Murad / Eur. J. Pure Appl. Math, 18 (3) (2025), 6633 22 of 24 of the problem (71-73) versus the numerical solution is depicted in Figures (6)a-b. Both graphs presented with N = 128, and the convergence achieved with tolerance 1e-06 at the 3rd iteration. (a) (b) Figure 5: Absolute error for the problem (68-70) as the mesh number (N) increases, (a) for a = 0.9 at t = 0.92, and (b) for a = 0.8 at t = 8.2. (a) (b) Figure 6: Exact solution (solid line) versus numerical solution (dashed line) for the problem (71-73), with µ = 0.1, α = β = 0.5, N = 128, (a) a = 0.9, and (b) a = 0.8. The reduction of absolute error is depicted in Figures (5)a-b as the mesh increases. For a = 0.9, the absolute error is of order O(1e − 04) and it is of order O(1e − 03) for a = 0.8. A similar behavior for the absolute error was observed for the previous example. 6. Conclusion This study has provided a detailed examination of the boundary-value problem (6-8), which is governed by the Lane-Emden fractional differential equation in terms of existence, uniqueness, and stability. The fractional derivative is of the Caputo-Hadamard type. The Banach and Sodavoski’s fixed-point theorems proved the existence and uniqueness of the solution. The stability of the problem of type Ulam-Hyers and Ulam-Hyers-Rassiass has been investigated. The existence, uniqueness, and stability of the solutions are demon- strated with an example in Section 4. The derived analytical solution to the problem is discretised on a graded mesh using the fractional rectangular, Lln,1 interpolation. The A. S. Hasan, S. A. Murad / Eur. J. Pure Appl. Math, 18 (3) (2025), 6633 23 of 24 resulting set of non-linear equations is solved using the Newton-Raphson method, which involves the Jacobian matrix to connect η(t) and θ(t, η(t)) in a non-linear way. The sta- bility and reliability of the proposed numerical scheme were studied through providing examples. References [1] I. Podlubny. Fractional differential equations: an introduction to fractional deriva- tives, fractional differential equations, to methods of their solution and some of their applications. elsevier, 1998. [2] Y. Zhou. Basic theory of fractional differential equations. World scientific, 2023. [3] B. Ahmad and R. P. Agarwal. Some new versions of fractional boundary value prob- lems with slit-strips conditions. Boundary Value Problems, 2014:1–12, 2014. [4] K. Diethelm and N. J. Ford. Analysis of fractional differential equations. Journal of Mathematical Analysis and Applications, 265(2):229–248, 2002. [5] Muhammad Farhan, Zahir Shah, Rashid Jan, and Saeed Islam. A fractional modeling approach of buruli ulcer in possum mammals. Physica Scripta, 98(6):065219, 2023. [6] A. S. Hasan. Numerical solution of the bagley-torvik equation using the integer-order derivatives expansion. Science Journal of University of Zakho, 6(2):64–69, 2018. [7] A. A. Kilbas. Theory and applications of fractional differential equations. North- Holland Mathematics Studies, 204, 2006. [8] Sh. A. Murad, H. J. Zekri, and S. Hadid. Existence and uniqueness theorem of fractional mixed volterra-fredholm integrodifferential equation with integral boundary conditions. International Journal of Differential Equations, 2011(1):304570, 2011. [9] H. J. Lane. On the theoretical temperature of the sun, under the hypothesis of a gaseous mass maintaining its volume by its internal heat, and depending on the laws of gases as known to terrestrial experiment. American Journal of Science, 2(148):57– 74, 1870. [10] R. Emden. Gaskugeln: Anwendungen der mechanischen Warmetheorie auf kosmolo- gische und meteorologische Probleme... BG Teubner, 1907. [11] S. Chandrasekhar. An introduction to the study of stellar structure, volume 2. Courier Corporation, 1957. [12] Dina Prialnik. An introduction to the theory of stellar structure and evolution. Cam- bridge University Press, 2009. [13] D. A. Frank-Kamenetskii. Diffusion and heat exchange in chemical kinetics, volume 2171. Princeton University Press, 2015. [14] R. W. Ibrahim. Existence of nonlinear lane-emden equation of fractional order. Miskolc Mathematical Notes, 13(1):39–52, 2012. [15] R. W. Ibrahim. Stability of a fractional differential equation. International Journal of Mathematical and Computational Sciences, 7(3):300–305, 2013. [16] A. K. Nasab, Z. P. Atabakan, A. I. Ismail, and R. W. Ibrahim. A numerical method for solving singular fractional lane–emden type equations. Journal of King Saud University-Science, 30(1):120–130, 2018. A. S. Hasan, S. A. Murad / Eur. J. Pure Appl. Math, 18 (3) (2025), 6633 24 of 24 [17] A. Taieb and Z. Dahmani. The hight order lane-emden fractional differential system: existence, uniqueness and ulam type stabilities. Kragujevac Journal of Mathematics, 40(2):238–259, 2016. [18] J. Dávila, L. Dupaigne, and J. Wei. On the fractional lane-emden equation. Trans- actions of the American Mathematical Society, 369(9):6087–6104, 2017. [19] M. E. Omaba. New analytical method of solution to a nonlinear singular fractional lane–emden type equation. AIMS Mathematics, 7(10):19539–19552, 2022. [20] Y. Gouari and Z. Dahmani. Stability of solutions for two classes of fractional dif- ferential equations of lane-emden type. Journal of Interdisciplinary Mathematics, 24(8):2087–2099, 2021. [21] R. O. Awonusika. Analytical solutions of a class of fractional lane–emden equation: A power series method. International Journal of Applied and Computational Mathe- matics, 8(4):155, 2022. [22] N. M. Dien. Solvability of nonlinear fractional lane–emden-type delay equations with time-singular coefficients. Rocky Mountain Journal of Mathematics, 54(3):855–868, 2024. [23] M. Gohar, C. Li, and Z. Li. Finite difference methods for caputo–hadamard fractional differential equations. Mediterranean Journal of Mathematics, 17(6):194, 2020. [24] Changpin Li, Zhiqiang Li, and Zhen Wang. Mathematical analysis and the local dis- continuous galerkin method for caputo–hadamard fractional partial differential equa- tion. Journal of Scientific Computing, 85:1–27, 2020. [25] C. W. H. Green, Y. Liu, and Y. Yan. Numerical methods for caputo–hadamard fractional differential equations with graded and non-uniform meshes. Mathematics, 9(21):2728, 2021. [26] F. Jarad, T. Abdeljawad, and D. Baleanu. Caputo-type modification of the hadamard fractional derivatives. Advances in Difference Equations, 2012:1–8, 2012. [27] Sh. Aljoudi, B. Ahmad, J. J. Nieto, and A. Alsaedi. A coupled system of hadamard type sequential fractional differential equations with coupled strip conditions. Chaos, Solitons & Fractals, 91:39–46, 2016. [28] S. Banach. Sur les opérations dans les ensembles abstraits et leur application aux équations intégrales. Fundamenta mathematicae, 3(1):133–181, 1922. [29] Ioan A Rus. Ulam stability of ordinary differential equations. Studia Universitatis Babeş-Bolyai, Mathematica, (4), 2009. [30] Saleh S Redhwan, Sadikali L Shaikh, Mohammed S Abdo, Wasfi Shatanawi, Ka- maleldin Abodayeh, Mohammed A Almalahi, and Tariq Aljaaidi. Investigating a generalized hilfer-type fractional differential equation with two-point and integral boundary conditions. AIMS Mathematics, 7(2):1856–1872, 2022. [31] Yi Yang and Jin Huang. Double fast algorithm for solving time-space fractional diffusion problems with spectral fractional laplacian. Applied Mathematics and Com- putation, 475:128715, 2024.