EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 4, Article Number 6388 ISSN 1307-5543 – ejpam.com Published by New York Business Global Qualitative Analysis and Simulation of Fractional Hybrid Boundary Value Problems in Orthogonal Cone Metric Spaces Dumitru Baleanu1,2, Mahammad Khuddush3, B.M.B. Krushna4, Sanket Tikare5,∗ 1 Department of Computer Science and Mathematics, Lebanese American University, Beirut, 1102 2801, Lebanon 2 Institute of Space Sciences-subsidiary of INFLPR, Magurele, Bucharest, 077125, Romania 3 Applied Nonlinear Science Lab(ANSL), Anand International College of Engineering, Jaipur 303012, India 4 Department of Mathematics, MVGR College of Engineering (Autonomous), Vizianagaram, 535 005, Andhra Pradesh, India 5 Department of Mathematics, Ramniranjan Jhunjhunwala College, Mumbai, 400 086, Maharashtra, India Abstract. This article aims to advance the qualitative analysis of fractional hybrid boundary value problems (FHBVPs) involving Riemann–Liouville fractional derivatives of order 1 < y ≤ 2, with a focus on establishing the existence, uniqueness, and stability of solutions in a novel math- ematical framework. By employing the extended Banach fixed point theorem within orthogonal cone metric spaces, we prove the existence and uniqueness of solutions for FHBVPs, generalizing prior results in standard metric spaces and Banach algebras [1, 2]. Additionally, we investigate the Hyers–Ulam stability to ensure solution robustness against perturbations, addressing com- mon methodological errors in prior studies [3]. Numerical simulations complement our theoretical findings, demonstrating the impact of fractional order and nonlinear terms on solution behavior. These results provide new insights into modeling complex dynamic systems with nonlocal and memory-dependent behaviors, applicable to fields such as viscoelasticity, fluid dynamics, and biological modeling. 2020 Mathematics Subject Classifications: 34A08, 47H10, 54H25 Key Words and Phrases: Existence and uniqueness, fractional derivative, Hyers–Ulam sta- bility, hybrid boundary value problem, orthogonal cone metric space ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.6388 Email addresses: dumitru.baleanu@lau.edu.lb (D. Baleanu), khuddush89@gmail.com (M. Khuddush), muraleebalu@yahoo.com (B.M.B. Krushna), sankettikare@rjcollege.edu.in (S. Tikare) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 2 of 31 1. Introduction Differential equations (DEqs) have long been a cornerstone of scientific and engineer- ing disciplines due to their ability to model dynamic systems and processes. However, classical DEqs often fail to capture phenomena with memory effects or nonlocal behavior. Fractional differential equations (FDEqs), which extend classical DEqs by incorporating derivatives of arbitrary order, address these limitations by modeling multiscale dynamics, making them particularly suitable for applications in viscoelasticity [4], fluid dynamics [5], and biological modeling [6–11]. The study of fractional calculus began in the 17th century with discussions by Leibniz and L’Hôpital on noninteger order derivatives, fol- lowed by formalization through the Riemann–Liouville fractional derivative in the 19th century [6]. The 20th century marked significant advancements, with FDEqs finding widespread applications in real-world problems. This evolution led to the development of fractional boundary value problems (FBVPs), which extend classical boundary value problems to model systems with nonlocal conditions and history-dependent behaviors [12, 13]. Concurrently, fixed point theory, particularly the Banach fixed point theorem, has been instrumental in establishing the existence and uniqueness of solutions for FB- VPs, with recent extensions to cone metric spaces [14] and orthogonal metric spaces [15] providing powerful tools for addressing complex nonlinear problems. Fractional hybrid boundary value problems (FHBVPs) further generalize FBVPs by incorporating quadratic perturbations of nonlinear DEqs, combining continuous and discrete dynamics. These equations are particularly valuable in applications such as biological systems, control theory, and economics, where sudden changes interact with continuous evolution [16–19]. However, prior studies on FHBVPs have notable limita- tions. For instance, Zhao et al. [1] investigated the existence and uniqueness of solutions for the Riemann–Liouville-based FHBVP: RLDα 0+ [ r(ω) f(ω, r(ω)) ] = g(ω, r(ω)) a.e. 0 < ω < ξ, 0 < α < 1, r(0) = 0. Their approach is constrained by its focus on lower-order fractional derivatives (0 < α < 1) and reliance on standard metric spaces with the Banach fixed point theorem, limiting its applicability to systems with higher-order derivatives or complex nonlocal boundary conditions. Similarly, Hilal and Kajouni [2] addressed the Caputo-type FHBVP: CDα 0+ [ r(ω) f(ω, r(ω)) ] = g(ω, r(ω)) a.e. 0 < ω < ξ, 0 < α < 1, a r(0) f(0, r(0)) + b r(ξ) f(ξ, r(ξ)) = c. Their framework, based on Banach algebra techniques and Lipschitz/Carathéodory con- ditions, is restricted by the use of Caputo derivative and specific algebraic structures, reducing its generality for broader classes of FHBVPs with higher-order derivatives or intricate boundary conditions. D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 3 of 31 To address these limitations, this study focuses on the qualitative analysis of the following fractional hybrid boundary value problem (FHBVP): RLDy 0+ [ r(ω) f(ω, r(ω)) ] + g(ω, r(ω)) = 0 a.e. 0 ≤ ω ≤ 1, (1) with boundary conditions (BCs) r(0) = 0 and r(1) = f ( 1, r(1) ) , (2) where 1 < y ≤ 2, RLDy 0+ is the Riemann–Liouville fractional derivative, g ∈ C([0, 1] × R,R), and f ∈ C([0, 1] × R,R \ {0}). We aim to establish the existence, uniqueness, and Hyers–Ulam stability of solutions using the extended Banach fixed point theorem in orthogonal cone metric spaces [20], which provides a more flexible geometric framework than standard metric spaces or Banach algebras. Recent advancements in metric spaces have significantly advanced fixed point theory. Traditional metric spaces have been extended to cone metric spaces by Huang et al. [14] and orthogonal metric spaces by Gordji et al. [15], enabling the analysis of complex nonlinear problems [21–24]. We know that the category of cone metric spaces and metric spaces are same and most fixed point results on cone metric spaces are not real generalizations. But, there are some valuable results (such results of this work) which researchers can work nowadays [25]. Our approach generalizes the results of Zhao et al. [1] and Hilal and Kajouni [2] by addressing higher-order fractional derivatives and leveraging the orthogonal cone metric space framework to accommodate intricate nonlinearities and nonlocal boundary conditions, offering new insights into the qualitative properties of FHBVPs across diverse scientific disciplines. The possible physical interpretations of (1)-(2): (i) Viscoelastic Material Deformation: • Meaning: The FHBVP can represent the deformation of a viscoelastic mate- rial (e.g., rubber, biological tissue) under a time-varying load. The fractional derivative RLDy 0+ captures the material’s memory-dependent stress relaxation, where y > 1 accounts for both elastic and viscous effects. The term r(ω) f(ω,r(ω)) might models a nonlinear stress-strain relationship adjusted by material prop- erties (f), and g(ω, r(ω)) could represent an external force or internal damp- ing. The boundary condition r(0) = 0 indicates no initial deformation, while r(1) = f(1, r(1)) suggests a terminal state dependent on the material’s non- linear response, possibly a fixed strain at the end of the interval. • Context: Applicable in engineering (e.g., designing shock absorbers) or biomechanics (e.g., modeling soft tissue). (ii) Anomalous Diffusion in Heterogeneous Media: • Meaning: This problem can describe anomalous diffusion processes, such as the spread of particles in porous media or fractals, where classical Fick’s law D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 4 of 31 fails. The fractional derivative reflects subdiffusion or superdiffusion due to memory effects over the interval [0, 1] (e.g., normalized time or space). The term r(ω) f(ω,r(ω)) might represents a concentration field scaled by a nonlinear permeability or reaction rate (f), and g(ω, r(ω)) could model source/sink terms (e.g., injection/extraction). The boundary conditions r(0) = 0 (no initial concentration) and r(1) = f(1, r(1)) (a nonlinear equilibrium at the boundary) suggest a system with controlled influx or outflow. • Context: Relevant in environmental science (e.g., groundwater flow) or ma- terial science (e.g., diffusion in nanocomposites). The novelty of this work lies in the following contributions: • Generalization to Higher-Order Fractional Derivatives: Unlike prior stud- ies such as Zhao et al. [1], which focus on FHBVPs with Riemann–Liouville deriva- tives of order 0 < α < 1, our work extends the analysis to higher-order derivatives (1 < y ≤ 2), enabling the modeling of more complex dynamic systems. • Orthogonal Cone Metric Space Framework: We apply the extended Banach fixed point theorem (FPT) in orthogonal cone metric spaces [20], generalizing the standard metric space and Banach algebra approaches of [1, 2]. This framework accommodates intricate nonlinearities and nonlocal boundary conditions, offering a more flexible geometric structure. • Hyers–Ulam Stability Analysis: Our study provides a rigorous Hyers–Ulam stability analysis for FHBVPs, correcting common misunderstandings in prior ap- proaches [3] by introducing a parameter to handle perturbed boundary conditions, ensuring robust solutions. • Numerical Validation: We complement our theoretical results with numerical simulations, using the trapezoidal rule to approximate solutions and providing graphical comparisons (e.g., Figures 1-4) to illustrate the impact of fractional order and function choices on solution behavior. The layout of this paper is as follows: Section 2 presents foundational definitions and key results required for later sections. In Section 3, we establish and prove several preliminary results that provide a basis for the main findings. Section 4 is dedicated to an in-depth examination of the existence and uniqueness of solutions for the FHBVP (1)–(2). Section 5 outlines the requirements for Hyers–Ulam stability concerning (1)– (2). Illustrative examples to support the obtained results are included in Section 6. In Section 7, we provide a conclusion of the work done in this paper and highlight possible avenues for future research. 2. Essential Preliminaries Definition 1. [6] For a continuous function Ψ : (0,∞) → R and an order y > 0, the D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 5 of 31 Riemann–Liouville fractional integral is described as follows: I y 0+ Ψ(ω) = 1 Γ(y) ∫ ω 0 (ω − ξ)y−1Ψ(ξ)dξ. Definition 2. [6] For a continuous function Ψ : (0,∞) → R and an order y > 0, the Riemann–Liouville fractional derivative is formulated as follows: Dy 0+ Ψ(ω) = 1 Γ(m− y) ( d dω )m ∫ ω 0 Ψ(ξ) (ω − ξ)y−m+1 dξ, m = [y] + 1. Remark 1. [6] The following composition relations are necessary for the present work: (i) Dy 0+ I y 0+ Ψ(ω) = Ψ(ω), y > 0, where Ψ(ω) ∈ L1(0,+∞). (ii) Dδ 0+I y 0+ Ψ(ω) = I y−δ 0+ Ψ(ω), y > δ > 0, where Ψ(ω) ∈ L1(0,+∞). Remark 2. [26] For α > −1, we have Dy 0+ ωα = Γ(α+ 1) Γ(α− y+ 1) ωα−y, which precisely yields Dy 0+ ωy−m = 0, m = 1, N, where N ≤ y ≤ N + 1 and N ∈ Z. Here 1, N = 1, 2, 3, . . . , N. Lemma 1. [6] Let y ∈ (m − 1, m] and m > 1. Then the general solution to Dy 0+ u(ω) = 0 is u(ω) = ∑m i=1 kiω y−i, where ki ∈ R, i = 1, m. Lemma 2. [6] Let y > 0. Then, for a given function u, we have I y 0+ Dy 0+ u(ω) = u(ω) + m∑ i=1 kiω y−i, where ki ∈ R, i = 1, m, m ≤ y ≤ m+ 1 and m ∈ Z. Definition 3. [14] Let E be a real Banach space and Q be a subset of E. Then Q is said to be a cone provided (i) Q is nonempty closed, and Q ̸= {0}; (ii) for l,m ∈ R, l,m ≥ 0, z, u ∈ Q, we have lz+mu ∈ Q; and (iii) z ∈ Q and −z ∈ Q imply z = 0. Definition 4. [14] Let A denote a nonempty collection of elements. Consider the func- tion γ : A×A → E which fulfills the following criteria: (c1) For all m,n ∈ A, it holds that 0 < γ(m,n), and γ(m,n) = 0 if and only if m = n. (c2) The equality γ(m,n) = γ(n,m) = 0 is satisfied for m,n ∈ A. D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 6 of 31 (c3) The inequality γ(m,n) ≤ γ(m, p) + γ(p, n) is true for m,n, p ∈ A. When the function γ satisfies these properties, it is referred to as a cone metric on A, and the pair (A, γ) is termed a cone metric space. Definition 5. [14] A cone C is termed normal provided there exists a positive constant K > 0 such that for every ℓ1, ℓ2 ∈ E, 0 ≤ ℓ1 ≤ ℓ2 implies ∥ℓ1∥ ≤ K∥ℓ2∥. The normality constant associated with C is the smallest positive value of K that satisfies this condition. Theorem 1. [14] Let (A, γ) denote a complete cone metric space, and let C be a normal cone characterized by its normality constant K. Suppose the function G : A → A fulfills the contractive property γ(Gu,Gv) ≤ cγ(u, v) for any u, v ∈ A, where c is a constant such that 0 ≤ c < 1. Given these assumptions, the mapping G possesses a unique fixed point within A. Additionally, for any element u ∈ A, the sequence defined by iterating G, denoted as {Gnu}, converges to this distinct fixed point. Definition 6. [15] Let A ̸= ∅ and ⊥⊆ A × A be a binary relation. The relation ⊥ is referred to as an orthogonal set (or simply an O-set) provided there exists an element α0 ∈ A such that for every β ∈ A, either β ⊥ α0 or α0 ⊥ β. We denote this O-set by (A,⊥). Definition 7. [15] Let (A,⊥) be an O-set. A sequence {αn}n∈N in A is termed an orthogonal sequence (or briefly an O-sequence) provided for all n ∈ N, either αn ⊥ αn+1 or αn ⊥ αn+1. Definition 8. [15] Let (A, γ,⊥) be an orthogonal metric space (where (A,⊥) is an O- set and (A, γ) is a metric space). The space A is said to be orthogonally complete (or briefly, O-complete) provided every Cauchy O-sequence converges. Remark 3. It is important to observe that while every complete metric space qualifies as O-complete, the reverse statement does not hold true (refer to [15]). Definition 9. [15] Let (A, γ,⊥) be an orthogonal metric space, and let 0 < c < 1. Then we have the following: 1. Orthogonal Contractive Function: A function G : A → A is called an orthog- onal contractive (or ⊥-contractive) function with Lipschitz constant c provided it satisfies γ(Gu,Gv) ≤ cγ(u, v), where u ⊥ v. 2. Orthogonal Preserving Function: A function G : A → A is termed an orthog- onal preserving (or ⊥-preserving) function provided G(u) ⊥ G(v) whenever u ⊥ v. 3. Orthogonally Continuous Function: A function G : A → A is classified as orthogonally continuous (or ⊥-continuous) at a point b ∈ A provided for any O- sequence {bn}n∈N in A, the convergence bn → b implies that G(bn) → G(b). Fur- thermore, G is regarded as ⊥-continuous on A provided it maintains ⊥-continuity at every b ∈ A. D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 7 of 31 Definition 10. [20] Let (A,⊥) be a nonempty orthogonal set. Assume the function γ : A×A → E satisfies the following conditions: (P1) 0 < γ(u, v) for all u, v ∈ A and γ(u, v) = 0 ⇐⇒ u = v; (P2) γ(u, v) = γ(v, u) = 0 for all u, v ∈ A; (P3) γ(u, v) ≤ γ(u,w) + γ(w, v) for all u, v, w ∈ A. Then γ is referred to as a cone metric on (A,⊥), and (A, γ,⊥) is classified as an orthogonal cone metric space. Definition 11. [20] Let (A, γ,⊥) be an orthogonal cone metric space. The space A is considered an orthogonally complete cone metric space provided every Cauchy O-sequence converges within A. Remark 4. It is important to note that every complete cone metric space qualifies as an O-complete space; however, the reverse is not necessarily valid (see [20]). Below we state the Banach FPT on an orthogonal cone metric space. Theorem 2. [20] Let (A, γ,⊥) be an O-complete metric space (which may not necessar- ily be complete). Consider a mapping G : A → A that satisfies the following conditions: 1. ⊥-continuity: The mapping G is continuous with respect to the orthogonal struc- ture present in the space. 2. ⊥-contraction: There exists a Lipschitz constant c such that for all u, v ∈ A, the inequality γ(G(u),G(v)) ≤ c · γ(u, v) holds, where c is constrained within the interval [0, 1). 3. ⊥-preserving: The mapping G maintains the orthogonality relation within the space. Given these criteria, the mapping G possesses a unique fixed point u∗ in the set A. Furthermore, G qualifies as a Picard operator, which implies that limr→∞ Gr(u) = u∗ for every u ∈ A. 3. Auxiliary Results In this section, we shall state and prove certain auxiliary results which serve as basis material for our main results. Lemma 3. The function r is a solution of the FHBVP (1)–(2) if and only if r is a solution to the integral equation r(ω) = f ( ω, r(ω) ) ωy−1 + f ( ω, r(ω) ) ∫ 1 0 0 ( ω, ξ ) g(ξ, r(ξ))dξ, (3) where 0(ω, ξ) = 1 Γ(y) { ωy−1(1− ξ)y−1 − (ω − ξ)y−1, ξ ≤ ω, ωy−1(1− ξ)y−1, ω ≤ ξ. (4) D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 8 of 31 Proof. Let r be a solution of the FHBVP (1)–(2). Then employing Lemma 2, we have r(ω) f ( ω, r(ω) ) = 2∑ i=1 aiω y−i − ∫ ω 0 (ω − ξ)y−1 Γ(y) g(ξ, r(ξ))dξ, (5) where a1, a2 are real constants. Applying the stipulations (2), it can be concluded that a2 = 0 as well as a1 = 1 + ∫ 1 0 (1− ξ)y−1 Γ(y) g(ξ, r(ξ))dξ. Now, putting a1 and a2 in (5), we get r(ω) f ( ω, r(ω) ) = ωy−1 + ∫ 1 0 ωy−1(1− ξ)y−1 Γ(y) g(ξ, r(ξ))dξ− ∫ ω 0 (ω − ξ)y−1 Γ(y) g(ξ, r(ξ))dξ. Therefore r(ω) = f ( ω, r(ω) ) ωy−1 + f ( ω, r(ω) ) ∫ 1 0 0 ( ω, ξ ) g(ξ, r(ξ))dξ. Assume, on the other hand, that r is a solution to (3). Then r(ω) f ( ω, r(ω) ) = ωy−1 + ∫ 1 0 0(ω, ξ)g(ξ, r(ξ))dξ = ωy−1 [ 1 + ∫ 1 0 (1− ξ)y−1 Γ(y) g(ξ, r(ξ))dξ ] − ∫ ω 0 (ω − ξ)y−1 Γ(y) g(ξ, r(ξ))dξ. Keeping in mind remarks 1 and 2, the operator RLD y 0+ is applied to each part of the aforementioned equation to arrive at RLDy 0+ [ r(ω) f(ω, r(ω)) ] = −g(ω, r(ω)). Given f(ω, r(ω)) ̸= 0 for ω ∈ [0, 1], it implies from (3) that r(0) = 0 and r(1) = f ( 1, r(1) ) . Consequently, the proof is now concluded. Lemma 4. The kernel 0(ω, ξ) possesses the following characteristics: (κ1) 0(ω, ξ) nonnegative and continuous on [0, 1]× [0, 1]. (κ2) 0(ω, ξ) ≤ 0(ξ, ξ) for ω, ξ ∈ [0, 1]. Proof. (κ1) From (4), we see that 0(ω, ξ) is continuous for ω, ξ ∈ [0, 1]. Next, for 0 ≤ ξ ≤ ω ≤ 1, we have 0(ω, ξ) = 1 Γ(y) [ ωy−1(1− ξ)y−1 − (ω − ξ)y−1 ] D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 9 of 31 = 1 Γ(y) [ ωy−1(1− ξ)y−1 − ωy−1 ( 1− ξ ω )y−1 ] ≥ 1 Γ(y) [ ωy−1(1− ξ)y−1 − ωy−1(1− ξ)y−1 ] = 0. Similarly, for 0 ≤ ω ≤ ξ ≤ 1, we find that 0(ω, ξ) ≥ 0. Thus, 0(ω, ξ) ≥ 0 for all ω, ξ ∈ [0, 1]. (κ2) For 0 ≤ ξ ≤ ω ≤ 1, we have ∂0(ω, ξ) ∂ω = 1 Γ(y) [ (y− 1)ωy−2(1− ξ)y−1 − (y− 1)(ω − ξ)y−2 ] = ωy−1 Γ(y− 1) [ (1− ξ)y−2 − ( 1− ξ ω )y−2 ] ≤ 0, 1 < y ≤ 2. This implies that 0(ω, ξ) is nonincreasing with respect to ω on [ξ, 1]. Hence, for 0 ≤ ξ ≤ ω ≤ 1, 0(ω, ξ) ≤ 0(ξ, ξ). Also, for 0 ≤ ω ≤ ξ ≤ 1, we have ∂0(ω, ξ) ∂ω = 1 Γ(y) [ (y− 1)ωy−2(1− ξ)y−1 ] ≥ 0, 1 < y ≤ 2, which implies that 0(ω, ξ) is nondecreasing with respect to ω on [0, ξ]. Hence, for 0 ≤ ω ≤ ξ ≤ 1, 0(ω, ξ) ≤ 0(ξ, ξ). Thus, we conclude that 0(ω, ξ) ≤ 0(ξ, ξ) for ω, ξ ∈ [0, 1]. This completes the proof. Next, we will proceed to outline and validate our next auxiliary finding, which acts as the central component of this paper. Theorem 3. Consider an orthogonal complete cone metric space (A, γ,⊥). Let there be a mapping F : A → A that maintains the orthogonality property and exhibits continuity with respect to this structure. Suppose F fulfills a specific contractive condition, referred to as γ(Fz, Fu) ≤ aγ(u, Fu) [ 1 + γ(z, Fz) ] 1 + γ(z, u) + bγ(z, u), applicable to any two elements z and u from A that are orthogonal, alongside constants a and b that fall within the range [0, 1), and satisfy the condition a+ b < 1. Under these stipulated conditions, it follows that F has exactly one fixed point in the set A. Proof. By the definition of orthogonality, there exists an element z0 ∈ A, where for every z ∈ A, either z is orthogonal to z0 or z0 is orthogonal to z. This leads us to D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 10 of 31 conclude that either z0 ⊥ F(z0) or F(z0) ⊥ z0. Now, we can define a sequence of elements as follows: z1 = F(z0), z2 = F(z1) = F2(z0), ... zn+1 = F(zn) = Fn(z0), n ∈ N. Then γ(zn, zn+1) = γ(Fzn−1, Fzn) ≤ aγ(zn, Fzn) [ 1 + γ(zn−1, Fzn−1) ] 1 + γ(zn−1, zn) + bγ(zn−1, zn) ≤ aγ(zn, zn+1) [ 1 + γ(zn−1, zn) ] 1 + γ(zn−1, zn) + bγ(zn−1, zn), which implies (1 − a)γ(zn, zn+1) ≤ bγ(zn−1, zn), n ∈ N. From this, taking r = b 1−a < 1, we can write that γ(zn, zn+1) ≤ rγ(zn−1, zn), ≤ · · · ≤ rnγ(z0, z1). For p, q ≥ 1, we have γ(zp, zp+q) ≤ γ(zp, zp+1) + γ(zp+1, zp+q) ≤ γ(zp, zp+1) + γ(zp+1, zp+2) + γ(zp+2, zp+q) ≤ γ(zp, zp+1) + γ(zp+1, zp+2) + γ(zp+2, zp+3) + · · · + γ(zp+q−3, zp+q−2) + γ(zp+q−2, zp+q−1) + γ(zp+q−1, zp+q) ≤ rpγ(z0, z1) + rp+1γ(z0, z1) + · · ·+ rp+q−1γ(z0, z1) ≤ rp 1− r γ(z0, z1). Thus, γ(zp, zp+q) ≤ rp 1−rγ(z0, z1). As p→ ∞, we deduce {zn} forms a Cauchy O-sequence. Since (A, γ) is a complete orthogonal cone metric space, we can find an element z⋆ ∈ A where zn converges to z⋆ as n→ ∞. Our next objective is to demonstrate that z⋆ serves as a fixed point of F. For this, γ(z⋆, Fz⋆) ≤ γ(z⋆, Fzn) + γ(Fzn, Fz ⋆) ≤ γ(z⋆, Fzn) + aγ(z⋆, Fz⋆) [ 1 + γ(zn, Fzn) ] 1 + γ(zn, z⋆) + bγ(zn, z ⋆) ≤ 1 + γ(zn, z ⋆) (1− a) + γ(zn, z⋆)− aγ(zn, zn+1) [ γ(z⋆, zn+1) + bγ(zn, z ⋆) ] D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 11 of 31 → 0 as n→ ∞. Hence Fz⋆ = z⋆, that is, z⋆ is a fixed point of F. Finally, we prove that z⋆ is a unique fixed point of F. If u⋆ is also a fixed point of F, i.e., Fu⋆ = u⋆. Then γ(u⋆, z⋆) = γ(Fu⋆, Fz⋆) ≤ aγ(z⋆, Fz⋆) [ 1 + γ(u⋆, Fu⋆) ] 1 + γ(u⋆, z⋆) + bγ(u⋆, z⋆) = bγ(u⋆, z⋆). Since b < 1, it follows that γ(u⋆, z⋆) = 0, which yields u⋆ = z⋆. This establishes the uniqueness of the fixed point of F and completes the proof. Remark 5. For a = 0, Theorem 3 reduces the Banach FPT on orthogonal cone metric space, i.e., Theorem 2. 4. Existence and Uniqueness of Solutions Consider the set A = {r ∈ C(I,R) : r(ω) ≥ 0 for almost every ω ∈ I}, where I := [0, 1]. Define the Banach space E = R and the subset Q = [0,∞) within E. We introduce a partial ordering ⪯ on E with respect to Q, where u ⪯ v if and only if v − u ∈ Q. Now, consider a functional γ : A×A → E given by γ(r1, r2) = sup ω∈I |r1(ω)− r2(ω)| for r1, r2 ∈ A. Theorem 4. Suppose the subsequent assertions hold true: (B1) A function ℘ ∈ L1(I,R) can be found satisfying |g(ω, r)− g(ω, u)| ≤ ℘(ω) ≤ |g(ω, 0)| for all ω ∈ I and r, u ∈ A. (B2) A constant 0 < b < 1 exists so that |f(ω, r)− f(ω, u)| ≤ b|r− u|[1 + |r− u|] 1 + 2M for ω ∈ I and r, u ∈ A, where M := ∫ 1 0 0(ω, ω)℘(ω) dω. Under these conditions, the FHBVP (1)–(2) has a unique solution. D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 12 of 31 Proof. Define an orthogonality relation ⊥ on the set A as follows z ⊥ u if and only if z(ω) · u(ω) ≥ 0 for almost all ω within the interval I. This relation demonstrates that (A, γ,⊥) satisfies the conditions for a cone metric space. Additionally, because each function z in A is continuous over a closed and bounded subset of Euclidean space, it attains a supremum in (A, γ,⊥). As a result, we conclude that (A, γ,⊥) is a complete space. Now, we introduce a function F : (A, γ,⊥) → (A, γ,⊥) defined by Fr(ω) = f(ω, r(ω))ωy−1 + f(ω, r(ω)) ∫ 1 0 0(ω, ξ)g(ξ, r(ξ))dξ, for all ω ∈ I. We observe that r ∈ A is a solution to FHBVP (1)–(2) if and only if r is a fixed point of F. Firstly, we prove that F is a self-mapping on z. To prove this, let ω ∈ I and r ∈ A. Then, Fr(ω) = f(ω, r(ω))ωy−1 + f(ω, r(ω)) ∫ 1 0 0(ω, ξ)g(ξ, r(ξ))dξ ≥ f(ω, r(ω)) ∫ 1 0 0(ω, ξ)g(ξ, r(ξ))dξ ≥ 0. (6) Therefore, we have F(z) ⊆ A. Next, we verify that the conditions of Theorem 2 are met. F is ⊥-preserving: Let r(ω) ⊥ u(ω) for all ω ∈ I. From (6), Fr(ω) ≥ 0 for all ω ∈ I, which implies that Fr ⊥ Fu, i.e., F is ⊥-preserving. F is ⊥-rational contraction: Let r, u ∈ A and r ⊥ u. Then, we have |Fr(ω)− Fu(ω)| ≤ |ω|y−1|f(ω, r(ω))− f(ω, u(ω))| + ∣∣∣∣f(ω, r(ω)) ∫ 1 0 0(ω, ξ)g(ξ, r(ξ))dξ − f(ω, u(ω)) ∫ 1 0 0(ω, ξ)g(ξ, u(ξ))dξ ∣∣∣∣ ≤ |ω|y−1|f(ω, r(ω))− f(ω, u(ω))|+ ∣∣∣∣f(ω, r(ω))∫ 1 0 0(ω, ξ)[g(ξ, r(ξ))− g(ξ, 0)]dξ − [ f(ω, u(ω)) ∫ 1 0 0(ω, ξ)[g(ξ, u(ξ))− g(ξ, 0)]dξ − f(ω, r(ω)) ∫ 1 0 0(ω, ξ)g(ξ, 0)dξ ] −f(ω, u(ω)) ∫ 1 0 0(ω, ξ)g(ξ, 0)dξ ∣∣∣∣ ≤ |ω|y−1|f(ω, r(ω))− f(ω, u(ω))| + ∣∣∣∣f(ω, r(ω))∫ 1 0 0(ω, ξ)[g(ξ, r(ξ))− g(ξ, 0)]dξ − f(ω, u(ω)) ∫ 1 0 0(ω, ξ)g(ξ, 0)dξ ∣∣∣∣ + ∣∣∣∣f(ω, u(ω))∫ 1 0 0(ω, ξ)[g(ξ, u(ξ))− g(ξ, 0)]dξ−f(ω, r(ω)) ∫ 1 0 0(ω, ξ)g(ξ, 0)dξ ∣∣∣∣ D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 13 of 31 ≤ |ω|y−1|f(ω, r(ω))− f(ω, u(ω))| + ∣∣∣∣f(ω, r(ω)) ∫ 1 0 0(ω, ξ)℘(ξ)dξ − f(ω, u(ω)) ∫ 1 0 0(ω, ξ)℘(ξ)dξ ∣∣∣∣ + ∣∣∣∣f(ω, r(ω)) ∫ 1 0 0(ω, ξ)℘(ξ)dξ − f(ω, u(ω)) ∫ 1 0 0(ω, ξ)℘(ξ)dξ ∣∣∣∣ ≤ |ω|y−1|f(ω, r(ω))− f(ω, u(ω))| + ∣∣∣∣f(ω, r(ω)) ∫ 1 0 0(ω, ξ)℘(ξ)dξ − f(ω, u(ω)) ∫ 1 0 0(ω, ξ)℘(ξ)dξ ∣∣∣∣ + ∣∣∣∣f(ω, r(ω)) ∫ 1 0 0(ω, ξ)℘(ξ)dξ − f(ω, u(ω)) ∫ 1 0 0(ω, ξ)℘(ξ)dξ ∣∣∣∣ ≤ |ω|y−1|f(ω, r(ω))− f(ω, u(ω))|+ 2|f(ω, r(ω))− f(ω, u(ω))| ∫ 1 0 0(ξ, ξ)℘(ξ)dξ ≤ |f(ω, r(ω))− f(ω, u(ω))|+ 2|f(ω, r(ω))− f(ω, u(ω))|M ≤ |f(ω, r(ω))− f(ω, u(ω))| (1 + 2M) ≤ b|r− u|[1 + |r− u|] 1 + 2M (1 + 2M) = b|r− u|[1 + |r− u|]. But, for any 0 < a < 1, we have |Fr− Fu|(1 + |r− u|)− a|r− Fr|(1 + |u− Fu|) ≤ |Fr− Fu| ≤ b|r− u| [ 1 + |r− u| ] which implies∣∣Fr(ω)− Fu(ω) ∣∣ ≤ a|r(ω)− Fr(ω)|(1 + |u(ω)− Fu(ω)|) 1 + |r(ω)− u(ω)| + b|r(ω)− u(ω)|. Taking supremum on both sides over ω, we get γ(Fr, Fu) ≤ aγ(r, Fr)[1 + γ(u, Fu)] 1 + γ(r, u) + bγ(r, u), (7) which shows that F is ⊥-rational contractive, since a+ b < 1. F is ⊥-continuous: Consider an O-sequence {rn} in z that converges to a point r ∈ A. Since F preserves the ⊥-orthogonality property, the sequence {F(rn)} also qualifies as an O-sequence. For any natural number n, setting r = rn, u = r, and a = 0 in Equation (7) yields ∣∣F(rn)− F(r) ∣∣ ≤ b|rn − r|. Taking the limit as n approaches infinity, it follows that F is ⊥-continuous. By applying Theorem 3, we conclude that r is the unique fixed point of F, which represents the solution for the FHBVP specified in (1)-(2). This completes the proof. D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 14 of 31 5. Hyers–Ulam Stability Analysis For some positive ε, consider the inequality∣∣∣∣RLDy 0+ [ r(ω) f(ω, r(ω)) ] + g(ω, r(ω)) ∣∣∣∣ ≤ ε (8) for ω ∈ [0, 1]. The fractional hybrid boundary value problem (FHBVP) (1)–(2) is susceptible to misunderstandings when applying Hyers–Ulam stability, as noted by Agarwal et al. [3]. They identify two main issues: (P1) treating the exact solution as fixed and independent of the approximate solution, and (P2) assuming the approximate solution satisfies the original boundary conditions, leading to invalid stability claims (see [3], Section 2.2.1). To address these, we modify the approach by introducing a parameter θ. Thus, the FHBVP is regarded as Hyers–Ulam stable provided for any r ∈ A satisfying Inequality (8), there exists a parameter θ = θ(r, ε) = r(1)− f(1, r(1)) and a corresponding solution u(ω, θ) ∈ A of the modified problem r(0) = 0 and r(1) = f(1, r(1)) + θ, (9) such that |r(ω)− u(ω, θ)| ≤ Kε (10) for some K > 0 independent of ε, aligning with the methodology proposed in [3] (Section 2.2.2). Remark 6. We say that r ∈ A is a solution of the Inequality (8) provided there exists a function ψ ∈ A, which depends upon r, such that |ψ(ω)| ≤ ε and RLDy 0+ [ r(ω) f(ω, r(ω)) ] + g(ω, r(ω)) = ψ(ω) for ω ∈ [0, 1]. (11) Note that r is not required to satisfy the original boundary condition r(1) = f(1, r(1)), avoiding the mistake (P2) highlighted in [3]. Lemma 5. Let r ∈ A be a solution of (8). Suppose sup ω∈[0,1] |f(ω, r(ω))| ≤ Q for some Q > 0, and N = ∫ 1 0 0(ω, ω) dω. Then the inequality∣∣∣∣r(ω)− f(ω, r(ω))ωy−1 − f(ω, r(ω)) ∫ 1 0 0(ω, ξ)g(ξ, r(ξ)) dξ ∣∣∣∣ ≤ NQε (12) holds for ω ∈ [0, 1]. D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 15 of 31 Proof. From Remark 6, there exists ψ(ω) with |ψ(ω)| ≤ ε such that (11) holds. The integral representation for a solution z of (11) with boundary condition z(0) = 0 and z(1) = f(1, z(1)) + θ (for some θ) is given by z(ω) = f(ω, z(ω))ωy−1+f(ω, z(ω)) ∫ 1 0 0(ω, ξ)g(ξ, z(ξ)) dξ−f(ω, z(ω)) ∫ 1 0 0(ω, ξ)ψ(ξ) dξ. Since r does not necessarily satisfy the original BC, we consider the perturbation. Sub- stituting ψ(ω) into the integral form and using the bound |ψ(ξ)| ≤ ε, we get∣∣∣∣r(ω)− f(ω, r(ω))ωy−1 − f(ω, r(ω)) ∫ 1 0 0(ω, ξ)g(ξ, r(ξ)) dξ ∣∣∣∣ ≤ |f(ω, r(ω))| ∫ 1 0 |0(ω, ξ)||ψ(ξ)| dξ. Using |0(ω, ξ)| ≤ 0(ξ, ξ) (by the definition of 0) and |ψ(ξ)| ≤ ε, we have∫ 1 0 |0(ω, ξ)||ψ(ξ)| dξ ≤ ε ∫ 1 0 0(ξ, ξ) dξ = Nε. Thus, ∣∣∣∣r(ω)− f(ω, r(ω))ωy−1 − f(ω, r(ω)) ∫ 1 0 0(ω, ξ)g(ξ, r(ξ)) dξ ∣∣∣∣ ≤ QNε, completing the proof. Theorem 5. Suppose (B1) and (B2) hold, where (B1) ensures the existence and unique- ness of u(ω, θ) for the modified FHBVP (1)–(9) for any θ, and (B2) provides Lipschitz conditions on f and g with constants Lf and Lg. Further, suppose NQ < 1, where Q = sup ω∈[0,1] |f(ω, r(ω))| and N = ∫ 1 0 0(ω, ω) dω. Then the FHBVP (1)–(2) is Hyers–Ulam stable. Proof. Assume r ∈ A fulfills the condition given in Inequality (8), and let u(ω, θ) ∈ A denote the unique solution to the modified FHBVP (1)–(9) with θ = r(1) − f(1, r(1)). In the view of Lemma 5, we have∣∣∣∣r(ω)− f(ω, r(ω))ωy−1 − f(ω, r(ω)) ∫ 1 0 0(ω, ξ)g(ξ, r(ξ)) dξ ∣∣∣∣ ≤ NQε. The exact solution u(ω, θ) satisfies u(ω, θ) = f(ω, u(ω, θ))ωy−1 + f(ω, u(ω, θ)) ∫ 1 0 0(ω, ξ)g(ξ, u(ξ, θ)) dξ. D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 16 of 31 Consider the difference |r(ω)− u(ω, θ)| ≤ ∣∣∣∣r(ω)− f(ω, r(ω))ωy−1 − f(ω, r(ω)) ∫ 1 0 0(ω, ξ)g(ξ, r(ξ)) dξ ∣∣∣∣ + ∣∣∣∣f(ω, r(ω))ωy−1 + f(ω, r(ω)) ∫ 1 0 0(ω, ξ)g(ξ, r(ξ)) dξ −f(ω, u(ω, θ))ωy−1 − f(ω, u(ω, θ)) ∫ 1 0 0(ω, ξ)g(ξ, u(ξ, θ)) dξ ∣∣∣∣ . Using the Lipschitz condition |f(ω, r) − f(ω, u)| ≤ LF|r − u| and |g(ξ, r) − g(ξ, u)| ≤ LG|r− u|, and bounding the integral term, we get∣∣∣∣f(ω, r(ω)) ∫ 1 0 0(ω, ξ)[g(ξ, r(ξ))− g(ξ, u(ξ, θ))] dξ ∣∣∣∣ ≤ QLGN|r− u|. Thus, |r(ω)− u(ω, θ)| ≤ NQε+ LFω y−1|r− u|+ QLGN|r− u|. Now, taking supremum over ω ∈ [0, 1], we get ∥r− u∥ ≤ NQε+M∥r− u∥, where M = supω(LFω y−1 + QLGN). Since NQ < 1, we have ∥r− u∥(1−M) ≤ NQε, which yields ∥r− u∥ ≤ NQ 1−M ε, where K = NQ 1−M > 0 if M < 1. Under (B2), ensure M < 1 (e.g., by bounding LF and LG). This completes the proof, correcting the fixed-solution error (P1) as per [3]. 6. Numerical Illustrations Example 1. Consider the FHBVP RLDy 0+ [ r(ω) f(ω, r(ω)) ] + g(ω, r(ω)) = 0, 0 < ω < 1, (13) r(0) = 0, r(1) = f(1, r(1)), (14) where y = 3 2 , g(ω, r) = 2(1 + ω) + (1 + ω) sin(r), and f(ω, r) = ω + 4 cos(r) 19 + 38M . D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 17 of 31 Let ℘(ω) = 2(1 + ω). To verify condition (B1), compute: g(ω, r)− g(ω, u) = [2(1 + ω) + (1 + ω) sin(r)]− [2(1 + ω) + (1 + ω) sin(u)] = (1 + ω)(sin(r)− sin(u)), so |g(ω, r)− g(ω, u)| = (1 + ω)| sin(r)− sin(u)|. Using the mean value theorem, we have | sin(r)− sin(u)| ≤ |r− u|, and thus |g(ω, r)− g(ω, u)| ≤ (1 + ω)|r− u| ≤ 2(1 + ω) = ℘(ω), since ω ∈ [0, 1]. Also, g(ω, 0) = 2(1 + ω) = ℘(ω), so (B1) holds. Next, compute the constants using Green’s function 0(ω, ξ) as 0(ω, ω) = ω 1 2 (1− ω) 1 2 Γ(32) , where Γ(32) = √ π 2 . Then N = ∫ 1 0 0(ω, ω)dω = 2√ π ∫ 1 0 ω 1 2 (1− ω) 1 2dω = 2√ π B ( 3 2 , 3 2 ) . Since B ( 3 2 , 3 2 ) = π 8 , N = 2√ π · π 8 = √ π 4 . Next, M = ∫ 1 0 0(ω, ω)℘(ω)dω = 2√ π ∫ 1 0 ω 1 2 (1− ω) 1 2 · 2(1 + ω)dω. But ∫ 1 0 ω 1 2 (1− ω) 1 2 (1 + ω)dω = ∫ 1 0 ω 1 2 (1− ω) 1 2dω + ∫ 1 0 ω 3 2 (1− ω) 1 2dω = B ( 3 2 , 3 2 ) +B ( 5 2 , 3 2 ) = π 8 + π 16 = 3π 16 . Hence M = 4√ π · 3π 16 = 12π 16 √ π = 3 √ π 4 . Let Q = 9 + 10M 19 + 38M = 9 + 10 · 3 √ π 4 19 + 38 · 3 √ π 4 ≈ 0.3207. D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 18 of 31 For condition (B2), compute |f(ω, r)− f(ω, u)| = 4 19 + 38M | cos(r)− cos(u)|, | cos(r)− cos(u)| ≤ |r− u|, |f(ω, r)− f(ω, u)| ≤ 4 19 + 38 · 3 √ π 4 |r− u| ≈ 4 69.5144 |r− u| ≈ 0.0575|r− u|, with b = 0.0575 < 1. Thus, by Theorem 4, the FHBVP (13)–(14) has a unique solution. Now, for Hyers–Ulam stability, we find the roots of − 4 19 z2 + 15 19 z− √ π 4 [ 18 + 15 √ π 10 + 15 √ π ] = 0, that is, −0.2105z2 + 0.7895z− 0.5400 = 0, which are z1 ≈ 0.9, z2 ≈ 2.85. Also, NQ ≈ √ π 4 · 0.3207 ≈ 0.14210648. Thus, by Theorem 5, the FHBVP is Hyers–Ulam stable. Numerical Approximation of the Solution: The integral equation for r is approx- imated numerically using the form (3). To compute the integral term∫ 1 0 0(ω, ξ)g(ξ, r(ξ))dξ, we employed the trapezoidal rule for numerical integration. The kernel function 0(ω, ξ) is computed at discrete values of ξ, and the integration is performed over the interval [0, 1] using the np.trapezoid() function from the Python numpy library. Creation of the Numerical Table: The functional value r(ω) is computed for 100 discrete values of ω in the interval [0, 1], starting from the initial approximation r(ω) = 0. The numerical solution is obtained for two different values of the fractional order y = 1.5 and y = 2. For each value of ω, the corresponding values of r(ω) are calculated iteratively by solving the integral equation (3), and the results are tabulated in Table 4, see Appendix. Graphical Representation of the Results: For better understand the behavior of r across different values of y, the computed values of r(ω) are plotted for both y = 1.5 and y = 2.0. The plots are generated using the Matplotlib library in Python, providing a visual comparison of the two scenarios. The first plot (Figure 1) shows the behavior of r for y = 1.5, while the second plot (Figure 2) compares the solutions for both y = 1.5 and y = 2.0. D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 19 of 31 Plot for y = 1.5: The graph of r(ω) for y = 1.5 is plotted in Figure 1, showing how the function evolves over the interval [0, 1]. The results indicate a smooth increase in r(ω) as ω approaches 1. Figure 1: Numerical values of r(ω) for y = 1.5. The plot shows how the solution r(ω) evolves smoothly as ω increases from 0 to 1, illustrating the behavior of the system for this fractional order. Comparative Plot for y = 1.5 and y = 2.0: The second plot compares the two solutions. For both values of y, the function r exhibits a similar overall trend; however, for y = 2.0, the values of r(ω) are slightly higher, indicating the effect of the increased fractional order on the solution behavior. This is given in Figure 2. D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 20 of 31 Figure 2: Comparison of the numerical values of r(ω) for y = 1.5 and y = 2.0. The plot shows that while both solutions follow a similar trend, the solution for y = 2.0 is slightly higher, indicating the impact of the increased fractional order on the behavior of the system. Example 2. Consider the FHBVP (1)–(2) with y = 1.5 and g(ω, r) = 2(1 + ω) + (1 + ω) sin(r). Let f1(ω, r) = ω + 4 cos(r) 19 + 38M and f1(ω, r) = eω + 4 sin(r) 19 + 38M . D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 21 of 31 Figure 3: Comparison of solutions of r(ω) using two different f functions: f1(ω, r) = ω + 4 cos(r) 19+38M and f2(ω, r) = eω + 4 sin(r) 19+38M , with y = 1.5 and g(ω, r) = 2(1 + ω) + (1 + ω) sin(r). The plot, in Figure 3, comparing the values of solutions r(ω) for the two f functions (f1 and f2) reveals several key observations. Both functions show a similar increasing trend in r(ω) as ω increases from 0 to 1, indicating a positive correlation. However, f1 (blue line) consistently exhibits a higher magnitude than f2 (orange line), highlighting the significant impact of the chosen f function on the resulting r(ω) values. While f1 produces a smoother curve, f2 introduces more fluctuations. This can be attributed to f1 including a term proportional to ω, while f2 relies on sin(r), which may have a lesser effect on magnitude. Overall, the plot illustrates how variations in the f function influence both the magnitude and shape of the solutions for r(ω). Example 3. Consider the FHBVP (1)-(2) with y = 1.5 and f(ω, r) = ω + 4 cos(r) 19 + 38M . Let g1(ω, r) = 2(1 + ω) + (1 + ω) sin(r) and g2(ω, r) = (1 − ω) cos(r). The plot, in Figure 4, compares the value of r(ω) using two different g functions (g1 and g2) while keeping the same f function. Both solutions show a similar increasing trend as ω moves from 0 to 1, indicating a positive correlation. However, g1 (blue line) yields consistently higher magnitudes than g2 (orange line). The curve for g1 is smoother with a steeper upward slope, while g2 introduces more fluctuations and a less pronounced increase. These differences arise from the terms in each g function, with g1 contributing positively D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 22 of 31 to magnitude through 2(1+ω) and (1+ω) sin(r), while g2 relies on (1−ω) cos(r), which can lead to lower magnitudes and more oscillations. Overall, the plot highlights how variations in the g function affect the solutions for r(ω). Figure 4: Comparison of solutions r(ω) using two different g functions: g1(ω, r) = 2(1 + ω) + (1 + ω) sin(r) and g2(ω, r) = (1− ω) cos(r), with y = 1.5 and f(ω, r) = 4 sin(r) 19+38M . Example 4. Consider the FHBVP (1)–(2) with fractional order y = 1.5 and f(ω, r) = 1 + ω + sin(r), g4(ω, r) = (1 + ω) + 0.1 r2. This choice introduces a quadratic nonlinearity with a small scaling factor to control growth. The conditions (B1)–(B2) are satisfied since |g4(ω, r)− g4(ω, u)| = 0.1 |r2 − u2| ≤ 0.2M |r− u| for some M > 0, and |f(ω, r)− f(ω, u)| = | sin(r)− sin(u)| ≤ |r− u|, ensuring Lipschitz continuity. Numerical solutions were obtained with the trapezoidal rule on a uniform grid of 200 points. Table 1 lists representative values of r(ω), and Figure 5 compares g4 against g1(ω, r) = 2(1 + ω) + (1 + ω) sin(r) from Example 1. We observe that g4 grows more moderately than g1, with values reaching about 3.22 at ω = 1, compared with 3.91 for g1. The quadratic term thus produces smoother profiles while still increasing steeply near ω = 1. D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 23 of 31 ω r(ω) (g1) r(ω) (g4) 0.00 0.0000 0.0000 0.20 1.1218 0.9875 0.40 2.0504 1.7512 0.60 2.7247 2.2987 0.80 3.3183 2.7703 1.00 3.9071 3.2210 Table 1: Example 4: Comparison of r(ω) for g1 and g4 with y = 1.5 and f(ω, r) = 1 + ω + sin(r). Figure 5: Example 4: Comparison of r(ω) for g1 (blue) and g4 (orange) under y = 1.5. The quadratic nonlinearity in g4 reduces overall growth compared with g1. Example 5. Consider (1)–(2) with y = 1.8 and f3(ω, r) = 1 + sin(r), g5(ω, r) = (1 + ω)r. Here, g5 is linear in r while f3 removes the direct ω dependence, yielding a smoother forcing term. Both (B1) and (B2) are satisfied, since |g5(ω, r)− g5(ω, u)| = (1 + ω)|r− u| ≤ 2|r− u|, and f3 is Lipschitz with constant Lf = 1. Numerical solutions are summarized in Table 2. Compared with Example 4, values of r(ω) under f3 are consistently smaller, reaching only 2.42 at ω = 1 versus 3.22 for f. This demonstrates that removing the ω term in f reduces growth. Figure 6 highlights the difference in magnitudes, with f3 producing smoother and more subdued solutions. D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 24 of 31 ω r(ω) (f, y = 1.5) r(ω) (f3, y = 1.8) 0.00 0.0000 0.0000 0.20 0.9875 0.3866 0.40 1.7512 0.9004 0.60 2.2987 1.4634 0.80 2.7703 1.9698 1.00 3.2210 2.4181 Table 2: Example 5: Comparison of r(ω) for f (y = 1.5) and f3 (y = 1.8) under g5(ω, r) = (1 + ω)r. Figure 6: Example 5: Comparison of r(ω) for f (y = 1.5, blue) and f3 (y = 1.8, orange). The removal of ω in f3 reduces solution magnitudes. Example 6. Finally, we compare our formulation with related works by Zhao et al. [1] and Hilal–Kajouni [2], both with α = 0.75. Zhao et al. [1] considered the Riemann– Liouville formulation RLDα 0+ [ r(ω) f(ω,r(ω)) ] = g(ω, r(ω)), r(0) = 0, with f(ω, r) = 1 + cos(r) and g(ω, r) = ω + sin(r). Hilal and Kajouni [2] studied the Caputo version, CDα 0+ [ r(ω) f(ω,r(ω)) ] = g(ω, r(ω)), a r(0) f(0,r(0)) + b r(1) f(1,r(1)) = c, with a = b = c = 1. Using the Caputo integral representation and enforcing the boundary condition at each iteration, we obtained the values in Table 3. D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 25 of 31 Compared with Example 4, both Zhao and Hilal solutions are smaller in magnitude, with Zhao reaching about 2.19 at ω = 1 versus 3.22 for Example 4. The Hilal profile is further constrained by the boundary condition, starting at 0 and attaining 1.95 at ω = 1. Figure 7 plots these curves on a logarithmic scale, clearly separating the profiles across orders of magnitude. The differences illustrate how fractional order (α = 0.75 vs. y = 1.5) and boundary conditions influence the growth and shape of solutions. ω r(ω) (Ex. 4) r(ω) (Zhao) r(ω) (Hilal) 0.00 0.0000 0.0000 0.0000 0.20 0.9875 2.1354 0.0000 0.40 1.7512 2.2120 2.0102 0.60 2.2987 2.2731 2.1198 0.80 2.7703 2.3258 2.2040 1.00 3.2210 2.1856 1.9523 Table 3: Example 6: Comparison of r(ω) for Example 4 (y = 1.5), Zhao et al. (α = 0.75), and Hilal–Kajouni (α = 0.75, Caputo + BC). Figure 7: Example 6: Log-scale comparison of r(ω) for Example 4 (y = 1.5) (blue), Zhao et al. (α = 0.75, orange), and Hilal–Kajouni (α = 0.75, green). Zhao and Hilal yield smaller magnitudes, while Hilal is further constrained near the boundary due to the algebraic condition. The numerical investigations highlight the sensitivity of fractional Hammerstein boundary value problems to both the choice of nonlinearities and the fractional order. In Example 4, the quadratic term in g4 moderated the growth of solutions compared with g1, while Example 5 showed that removing the explicit ω-dependence in f3 reduced D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 26 of 31 the overall amplitude of solutions. Example 6 demonstrated that lowering the fractional order from y = 1.5 to α = 0.75 significantly decreased the magnitude of solutions, with Hilal’s Caputo formulation further constrained by the algebraic boundary condi- tion. Tables 1–3 and figures 5–7 together confirm that nonlinear structure, fractional order, and boundary conditions jointly govern the qualitative and quantitative behavior of solutions. 7. Conclusions The primary aim of this study is to develop a robust mathematical framework for an- alyzing fractional hybrid boundary value problems (FHBVPs) with Riemann–Liouville fractional derivatives of order 1 < y ≤ 2, focusing on establishing the existence, unique- ness, and stability of solutions. By leveraging the extended Banach fixed point theorem in orthogonal cone metric spaces, we successfully proved the existence and uniqueness of solutions for FHBVPs, extending prior results in standard metric spaces and Banach algebras [1, 2]. Our investigation of Hyers–Ulam stability, incorporating a parameter to address perturbed boundary conditions, ensures the robustness of solutions and corrects methodological issues noted in the literature [3]. Numerical simulations, implemented via the trapezoidal rule, validated our theoretical findings by illustrating the influence of fractional order and nonlinear terms on solution behavior, as demonstrated in figures 1–4. These results advance the theoretical understanding of FHBVPs and offer practical tools for modeling complex systems with nonlocal and memory-dependent dynamics in fields such as engineering, physics, and biology. Future research directions include the following: • Higher-Order and Multivariable Systems: Extend the analysis to FHBVPs with fractional orders beyond 1 < y ≤ 2 or multivariable hybrid systems to capture richer dynamics in applications like control theory and biological networks. • Diverse Boundary Conditions: Investigate the impact of alternative boundary conditions, such as multi-point or integral conditions, on the existence, uniqueness, and stability of solutions to broaden the applicability of the framework. • Advanced Numerical Methods: Develop more sophisticated numerical tech- niques, such as adaptive algorithms or machine learning-based approaches, to en- hance the accuracy and scalability of simulations for large-scale FHBVPs. • Real-World Applications: Apply the proposed framework to specific problems in viscoelasticity, anomalous diffusion, or biological systems, validating the model with experimental data to bridge theoretical and practical domains. • Generalized Metric Spaces: Explore other generalized metric spaces, such as partial metric spaces or fuzzy metric spaces, to further extend the fixed point techniques for FHBVPs and related problems. These avenues promise to deepen the understanding of fractional hybrid systems and enhance their applicability across diverse scientific disciplines. D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 27 of 31 Acknowledgements We would like to thank the reviewers for their valuable comments and constructive feedback, which have significantly contributed to the improvement of this manuscript. The authors would like to express their sincere gratitude to the respective institu- tions and collaborators for their support in this research. Specifically, M. Khuddush is thankful to the Applied Nonlinear Science Lab (ANSL) for providing the necessary research facilities and environment at AICE, Jaipur, India. B. M. B. Krushna is thank- ful to MVGR College of Engineering, Vizianagaram, India, for the support during the preparation of this paper. References [1] Y Zhao, S Sun, Z Han, and Q Li. Theory of fractional hybrid differential equations. Comput. Math. Appl., 62(3):1312–1324, 2011. [2] K Hilal and A Kajouni. Boundary value problems for hybrid differential equations with fractional order. Adv. Difference Equ., 2015(183):1–19, 2015. [3] RP Agarwal, S Hristova, and D O’Regan. Ulam stability for boundary value prob- lems of differential equations—main misunderstandings and how to avoid them. Mathematics, 12:1626, 2024. [4] TM Atanackovic, S Pilipovic, and B Stankovic D Zorica. Fractional calculus with applications in mechanics: Vibrations and diffusion processes. John Wiley & Sons, NJ, USA, 2014. [5] WG Glöckle and TF Nonnenmacher. A fractional calculus approach to self-similar protein dynamics. Biophys J., 68(1):46–53, 1995. [6] AA Kilbas, HM Srivastava, and JJ Trujillo. Theory and applications of fractional differential equations. Elsevier BV, Amsterdam, The Netherlands, 2006. [7] RL Magin. Fractional calculus models of complex dynamics in biological tissues. Comput. Math. Appl., 59(5):1586–1593, 2010. [8] SR Manam. Multiple integral equations arising in the theory of water waves. Appl. Math. Lett., 24(8):1369–1373, 2011. [9] KS Miller and B Ross. An Introduction to the Fractional Calculus and Fractional Differential Equations. John Wiley & Sons, New York, USA, 1993. [10] H Mohammadi, MKA Kaabar, J Alzabut, AGM Selvam, and S Rezapour. A com- plete model of crimean-congo hemorrhagic fever (cchf) transmission cycle with non- local fractional derivative. J. Funct. Spaces, 2021(1):1273405, 2021. [11] KR Prasad, BMB Krushna, VVRRB Raju, and Y Narasimhulu. Existence of posi- tive solutions for systems of fractional order boundary value problems with riemann– liouville derivative. Nonlinear Stud., 24(3):619–629, 2017. [12] M Feng, X Zhang, and W Ge. New existence results for higher-order nonlinear fractional differential equation with integral boundary conditions. Bound. Value Probl., 720702:1–20, 2011. [13] FJ Torres. Existence of a positive solution for a boundary value problem of a D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 28 of 31 nonlinear fractional differential equation. Bull. Iran. Math. Soc., 39(2):307–323, 2013. [14] LG Huang and X Zhang. Cone metric spaces and fixed point theorems of contractive mappings. J. Math. Anal. Appl., 332(2):1468–1476, 2007. [15] ME Gordji, M Ramezani, M De La Sen, and YJ Cho. On orthogonal sets and Banach fixed point theorem. Fixed Point Theory, 18(2):569–578, 2017. [16] B Ahmad, SK Ntouyas, and J Tariboon. A nonlocal hybrid boundary value problem of Caputo fractional integro-differential equations. Acta Math. Sci., 36:1631–1640, 2016. [17] BC Dhage and V Lakshmikantham. Basic results on hybrid differential equations. Nonlinear Anal. Hybrid Syst., 4(3):414–424, 2010. [18] AEM Herzallah and D Baleanu. On fractional order hybrid differential equations. Abstr. Appl. Anal., 2014(1):389386, 2014. [19] Z Ullah, A Ali, RA Khan, and M Iqbal. Existence results to a class of hybrid fractional differential equations. Matriks Sains Mat., 2(1):13–17, 2018. [20] ZEDD Olia, ME Gordji, and DE Bagha. Banach fixed point theorem on orthogonal cone metric spaces. Facta Univ., Ser. Math. Inf., 35(5):1239–1250, 2020. [21] LB Ćirić. A generalization of Banach’s contraction principle. Proc. Amer. Math. Soc., 45(2):267–273, 1974. [22] BK Lahiri, P Das, and LK Dey. Cantor’s theorem in 2-metric spaces and its appli- cations to fixed point problems. Taiwanese J. Math., 15(1):337–352, 2011. [23] KR Prasad, M Khuddush, and D Leela. Existence, uniqueness and hyers–ulam stability of a fractional order iterative two-point boundary value problems. Afr. Mat., 32:1227–1237, 2021. [24] T Senapati, LK Dey, and D Dolicanin-Dekic. Extension of ciric and wardowski type fixed point theorems in d-generalized metric spaces. Fixed Point Theory Appl., 2016(33):1–14, 2016. [25] S. Rezapour. Lights and Shadows on Generalizations in Fixed Point Theory. CRC Press, Taylor and Francis, 2025. [26] Z Bai and H Lü. Positive solutions for a boundary value problem of nonlinear fractional differential equations. J. Math. Anal. Appl., 311(2):495–505, 2005. Appendix Table 4, compares the values of r(ω) for y = 1.5 and y = 2.0 given in Example 1. Table 5, compares the solutions r(ω) for both functions f1 and f2 given in Example 2. Table 6, compares the solutions r(ω) for both functions g1 and g2 given in Example 3. D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 29 of 31 ω r(ω) (y = 1.5) r(ω) (y = 2.0) ω r(ω) (y = 1.5) r(ω) (y = 2.0) ω r(ω) (y = 1.5) r(ω) (y = 2.0) 0.00 0.0000 0.0000 0.34 0.5919 0.2686 0.68 1.0949 0.7474 0.01 0.0210 0.0016 0.35 0.6109 0.2817 0.69 1.1037 0.7614 0.02 0.0340 0.0036 0.36 0.6298 0.2949 0.70 1.1120 0.7754 0.03 0.0468 0.0061 0.37 0.6486 0.3084 0.71 1.1198 0.7891 0.04 0.0600 0.0091 0.38 0.6673 0.3220 0.72 1.1270 0.8027 0.05 0.0736 0.0125 0.39 0.6858 0.3357 0.73 1.1337 0.8161 0.06 0.0877 0.0163 0.40 0.7042 0.3497 0.74 1.1397 0.8292 0.07 0.1022 0.0205 0.41 0.7225 0.3637 0.75 1.1452 0.8422 0.08 0.1173 0.0252 0.42 0.7405 0.3779 0.76 1.1502 0.8549 0.09 0.1327 0.0302 0.43 0.7584 0.3923 0.77 1.1545 0.8674 0.10 0.1486 0.0357 0.44 0.7760 0.4067 0.78 1.1581 0.8796 0.11 0.1648 0.0416 0.45 0.7934 0.4213 0.79 1.1612 0.8916 0.12 0.1815 0.0478 0.46 0.8106 0.4359 0.80 1.1636 0.9033 0.13 0.1984 0.0545 0.47 0.8276 0.4507 0.81 1.1653 0.9146 0.14 0.2156 0.0615 0.48 0.8443 0.4655 0.82 1.1664 0.9257 0.15 0.2331 0.0689 0.49 0.8607 0.4804 0.83 1.1668 0.9365 0.16 0.2509 0.0766 0.50 0.8768 0.4953 0.84 1.1665 0.9469 0.17 0.2689 0.0847 0.51 0.8926 0.5103 0.85 1.1655 0.9570 0.18 0.2871 0.0932 0.52 0.9081 0.5253 0.86 1.1638 0.9667 0.19 0.3055 0.1019 0.53 0.9233 0.5404 0.87 1.1613 0.9761 0.20 0.3241 0.1110 0.54 0.9382 0.5554 0.88 1.1581 0.9850 0.21 0.3428 0.1205 0.55 0.9527 0.5705 0.89 1.1542 0.9936 0.22 0.3617 0.1302 0.56 0.9668 0.5855 0.90 1.1495 1.0018 0.23 0.3807 0.1403 0.57 0.9806 0.6006 0.91 1.1440 1.0095 0.24 0.3997 0.1506 0.58 0.9940 0.6156 0.92 1.1377 1.0168 0.25 0.4189 0.1613 0.59 1.0070 0.6305 0.93 1.1306 1.0236 0.26 0.4381 0.1722 0.60 1.0195 0.6454 0.94 1.1227 1.0300 0.27 0.4573 0.1834 0.61 1.0317 0.6603 0.95 1.1140 1.0358 0.28 0.4766 0.1948 0.62 1.0434 0.6751 0.96 1.1044 1.0412 0.29 0.4959 0.2065 0.63 1.0546 0.6897 0.97 1.0940 1.0461 0.30 0.5151 0.2185 0.64 1.0654 0.7043 0.98 1.0827 1.0505 0.31 0.5344 0.2307 0.65 1.0757 0.7188 0.99 1.0706 1.0543 0.32 0.5536 0.2431 0.66 1.0856 0.7331 1.00 1.0575 1.0575 0.33 0.5728 0.2558 0.67 1.0949 0.7474 Table 4: Comparison of r(ω) for y = 1.5 and y = 2.0 D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 30 of 31 ω (f1) (f2) ω (f1) (f2) ω (f1) (f2) 0.00 0.0000 0.0000 0.34 2.1659 0.1741 0.68 3.0195 0.6829 0.01 0.3320 0.0000 0.35 2.2018 0.1857 0.69 3.0321 0.6995 0.02 0.4715 0.0002 0.36 2.2372 0.1977 0.70 3.0438 0.7159 0.03 0.5800 0.0005 0.37 2.2720 0.2101 0.71 3.0544 0.7322 0.04 0.6727 0.0010 0.38 2.3063 0.2227 0.72 3.0639 0.7482 0.05 0.7554 0.0017 0.39 2.3400 0.2357 0.73 3.0724 0.7640 0.06 0.8310 0.0027 0.40 2.3731 0.2491 0.74 3.0798 0.7796 0.07 0.9014 0.0040 0.41 2.4057 0.2627 0.75 3.0860 0.7949 0.08 0.9677 0.0055 0.42 2.4377 0.2766 0.76 3.0911 0.8098 0.09 1.0306 0.0074 0.43 2.4691 0.2909 0.77 3.0949 0.8244 0.10 1.0907 0.0096 0.44 2.5000 0.3054 0.78 3.0976 0.8387 0.11 1.1485 0.0121 0.45 2.5302 0.3201 0.79 3.0990 0.8526 0.12 1.2043 0.0149 0.46 2.5599 0.3352 0.80 3.0991 0.8661 0.13 1.2583 0.0181 0.47 2.5889 0.3504 0.81 3.0978 0.8791 0.14 1.3107 0.0217 0.48 2.6173 0.3659 0.82 3.0953 0.8916 0.15 1.3618 0.0256 0.49 2.6451 0.3816 0.83 3.0913 0.9037 0.16 1.4116 0.0299 0.50 2.6722 0.3975 0.84 3.0859 0.9152 0.17 1.4602 0.0346 0.51 2.6987 0.4136 0.85 3.0791 0.9261 0.18 1.5078 0.0397 0.52 2.7245 0.4299 0.86 3.0708 0.9364 0.19 1.5544 0.0451 0.53 2.7496 0.4463 0.87 3.0609 0.9461 0.20 1.6000 0.0510 0.54 2.7740 0.4629 0.88 3.0495 0.9552 0.21 1.6449 0.0572 0.55 2.7977 0.4796 0.89 3.0365 0.9636 0.22 1.6889 0.0638 0.56 2.8207 0.4964 0.90 3.0218 0.9712 0.23 1.7322 0.0709 0.57 2.8429 0.5133 0.91 3.0055 0.9781 0.24 1.7748 0.0783 0.58 2.8643 0.5302 0.92 2.9875 0.9841 0.25 1.8167 0.0861 0.59 2.8850 0.5473 0.93 2.9676 0.9894 0.26 1.8579 0.0944 0.60 2.9048 0.5643 0.94 2.9460 0.9938 0.27 1.8985 0.1030 0.61 2.9238 0.5814 0.95 2.9225 0.9973 0.28 1.9384 0.1120 0.62 2.9420 0.5984 0.96 2.8971 0.9998 0.29 1.9778 0.1214 0.63 2.9593 0.6155 0.97 2.8698 1.0014 0.30 2.0165 0.1312 0.64 2.9757 0.6325 0.98 2.8405 1.0020 0.31 2.0547 0.1414 0.65 2.9913 0.6494 0.99 2.8092 1.0015 0.32 2.0923 0.1519 0.66 3.0058 0.6662 1.00 2.7758 1.0000 0.33 2.1294 0.1628 0.67 3.0058 0.6662 Table 5: Comparison of r(ω) for f1(ω, r) = ω + 4 cos(r) 19+38M and f2(ω, r) = eω + 4 sin(r) 19+38M D. Baleanu et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6388 31 of 31 ω (g1) (g2) ω (g1) (g2) ω (g1) (g2) 0.00 0.0000 0.0000 0.34 0.0087 0.0041 0.68 0.0088 0.0051 0.01 0.0018 0.0008 0.35 0.0088 0.0042 0.69 0.0087 0.0051 0.02 0.0026 0.0012 0.36 0.0088 0.0042 0.70 0.0086 0.0051 0.03 0.0032 0.0014 0.37 0.0088 0.0043 0.71 0.0086 0.0052 0.04 0.0036 0.0016 0.38 0.0089 0.0043 0.72 0.0085 0.0052 0.05 0.0040 0.0018 0.39 0.0089 0.0043 0.73 0.0085 0.0052 0.06 0.0044 0.0020 0.40 0.0090 0.0044 0.74 0.0084 0.0052 0.07 0.0047 0.0021 0.41 0.0091 0.0044 0.75 0.0083 0.0052 0.08 0.0050 0.0023 0.42 0.0091 0.0044 0.76 0.0083 0.0053 0.09 0.0053 0.0024 0.43 0.0091 0.0045 0.77 0.0082 0.0053 0.10 0.0056 0.0025 0.44 0.0091 0.0045 0.78 0.0081 0.0053 0.11 0.0058 0.0026 0.45 0.0092 0.0045 0.79 0.0080 0.0053 0.12 0.0060 0.0027 0.46 0.0092 0.0046 0.80 0.0079 0.0053 0.13 0.0062 0.0028 0.47 0.0092 0.0046 0.81 0.0079 0.0054 0.14 0.0064 0.0029 0.48 0.0092 0.0046 0.82 0.0078 0.0054 0.15 0.0066 0.0030 0.49 0.0092 0.0047 0.83 0.0077 0.0054 0.16 0.0068 0.0031 0.50 0.0092 0.0047 0.84 0.0076 0.0054 0.17 0.0070 0.0032 0.51 0.0092 0.0047 0.85 0.0075 0.0054 0.18 0.0071 0.0032 0.52 0.0092 0.0047 0.86 0.0074 0.0055 0.19 0.0073 0.0033 0.53 0.0092 0.0047 0.87 0.0073 0.0055 0.20 0.0074 0.0034 0.54 0.0092 0.0048 0.88 0.0072 0.0055 0.21 0.0075 0.0035 0.55 0.0092 0.0048 0.89 0.0071 0.0055 0.22 0.0077 0.0035 0.56 0.0092 0.0048 0.90 0.0070 0.0055 0.23 0.0078 0.0036 0.57 0.0092 0.0048 0.91 0.0069 0.0056 0.24 0.0079 0.0036 0.58 0.0091 0.0048 0.92 0.0068 0.0056 0.25 0.0080 0.0037 0.59 0.0091 0.0049 0.93 0.0066 0.0056 0.26 0.0081 0.0037 0.60 0.0091 0.0049 0.94 0.0065 0.0056 0.27 0.0082 0.0038 0.61 0.0090 0.0049 0.95 0.0064 0.0056 0.28 0.0083 0.0039 0.62 0.0090 0.0049 0.96 0.0063 0.0056 0.29 0.0084 0.0039 0.63 0.0090 0.0049 0.97 0.0062 0.0056 0.30 0.0085 0.0040 0.64 0.0090 0.0049 0.98 0.0061 0.0056 0.31 0.0086 0.0040 0.65 0.0089 0.0049 0.99 0.0060 0.0057 0.32 0.0087 0.0041 0.66 0.0089 0.0049 1.00 0.0059 0.0057 0.33 0.0087 0.0041 0.67 0.0088 0.0050 Table 6: Values of r(ω) for ω = 0 to ω = 1.00 Introduction Essential Preliminaries Auxiliary Results Existence and Uniqueness of Solutions Hyers–Ulam Stability Analysis Numerical Illustrations Conclusions