diff --git a/Project.toml b/Project.toml
index 8e603699ba..1da7238259 100644
--- a/Project.toml
+++ b/Project.toml
@@ -44,6 +44,7 @@ NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
@@ -141,6 +142,7 @@ OrdinaryDiffEq = "6.82.0"
 OrdinaryDiffEqCore = "1.15.0"
 OrdinaryDiffEqDefault = "1.2"
 OrdinaryDiffEqNonlinearSolve = "1.5.0"
+Plots = "1.40.15"
 PrecompileTools = "1"
 Pyomo = "0.1.0"
 REPL = "1"
diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index 1dbf9d37dd..dee6a22fcb 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -266,8 +266,9 @@ export isinput, isoutput, getbounds, hasbounds, getguess, hasguess, isdisturbanc
        tunable_parameters, isirreducible, getdescription, hasdescription,
        hasunit, getunit, hasconnect, getconnect,
        hasmisc, getmisc, state_priority
-export liouville_transform, change_independent_variable, substitute_component,
-       add_accumulations, noise_to_brownians, Girsanov_transform, change_of_variables
+export liouville_transform, change_independent_variable, substitute_component, add_accumulations,
+       noise_to_brownians, Girsanov_transform, change_of_variables,
+       fractional_to_ordinary, linear_fractional_to_ordinary
 export PDESystem
 export Differential, expand_derivatives, @derivatives
 export Equation, ConstrainedEquation
diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index b8268e884e..c077a7ffe2 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -185,6 +185,232 @@ function change_of_variables(
     return new_sys
 end
 
+"""
+Generates the system of ODEs to find solution to FDEs.
+
+Example:
+
+```julia
+@independent_variables t
+@variables x(t)
+D = Differential(t)
+tspan = (0., 1.)
+
+α = 0.5
+eqs = (9*gamma(1 + α)/4) - (3*t^(4 - α/2)*gamma(5 + α/2)/gamma(5 - α/2))
+eqs += (gamma(9)*t^(8 - α)/gamma(9 - α)) + (3/2*t^(α/2)-t^4)^3 - x^(3/2)
+sys = fractional_to_ordinary(eqs, x, α, 10^-7, 1)
+
+prob = ODEProblem(sys, [], tspan)
+sol = solve(prob, radau5(), abstol = 1e-10, reltol = 1e-10)
+
+time = 0
+while(time <= 1)
+    @test isapprox((3/2*time^(α/2) - time^4)^2, sol(time, idxs=x), atol=1e-3)
+    time += 0.1
+end
+```
+"""
+function fractional_to_ordinary(
+        eqs, variables, alphas, epsilon, T;
+        initials = 0, additional_eqs = [], iv = only(@independent_variables t), matrix=false
+)
+    D = Differential(iv)
+    i = 0
+    all_eqs = Equation[]
+    all_def = Pair[]
+    
+    function fto_helper(sub_eq, sub_var, α; initial=0)
+        alpha_0 = α
+
+        if (α > 1)
+            coeff = 1/(α - 1)
+            m = 2
+            while (α - m > 0)
+                coeff /= α - m
+                m += 1
+            end
+            alpha_0 = α - m + 1
+        end
+
+        δ = (gamma(alpha_0+1) * epsilon)^(1/alpha_0)
+        a = pi/2*(1-(1-alpha_0)/((2-alpha_0) * log(epsilon^-1)))
+        h = 2*pi*a / log(1 + (2/epsilon * (cos(a))^(alpha_0 - 1)))
+
+        x_sub = (gamma(2-alpha_0) * epsilon)^(1/(1-alpha_0))
+        x_sup = -log(gamma(1-alpha_0) * epsilon)
+        M = floor(Int, log(x_sub / T) / h)
+        N = ceil(Int, log(x_sup / δ) / h)
+
+        function c_i(index)
+            h * sin(pi * alpha_0) / pi * exp((1-alpha_0)*h*index)
+        end
+
+        function γ_i(index)
+            exp(h * index)
+        end
+
+        new_eqs = Equation[]
+        def = Pair[]
+
+        if matrix
+            new_z = Symbol(:ʐ, :_, i)
+            i += 1
+            γs = diagm([γ_i(index) for index in M:N-1])
+            cs = [c_i(index) for index in M:N-1]
+
+            if (α < 1)
+                new_z = only(@variables $new_z(iv)[1:N-M])
+                new_eq = D(new_z) ~ -γs*new_z .+ sub_eq
+                rhs = dot(cs, new_z) + initial
+                push!(def, new_z=>zeros(N-M))
+            else
+                new_z = only(@variables $new_z(iv)[1:N-M, 1:m])
+                new_eq = D(new_z) ~ -γs*new_z + hcat(fill(sub_eq, N-M, 1), collect(new_z[:, 1:m-1]*diagm(1:m-1)))
+                rhs = coeff*sum(cs[i]*new_z[i, m] for i in 1:N-M)
+                for (index, value) in enumerate(initial)
+                    rhs += value * iv^(index - 1) / gamma(index)
+                end
+                push!(def, new_z=>zeros(N-M, m))
+            end
+            push!(new_eqs, new_eq)
+        else
+            if (α < 1)  
+                rhs = initial
+                for index in range(M, N-1; step=1)
+                    new_z = Symbol(:ʐ, :_, i)
+                    i += 1
+                    new_z = ModelingToolkit.unwrap(only(@variables $new_z(iv)))
+                    new_eq = D(new_z) ~ sub_eq - γ_i(index)*new_z
+                    push!(new_eqs, new_eq)
+                    push!(def, new_z=>0)
+                    rhs += c_i(index)*new_z
+                end
+            else
+                rhs = 0
+                for (index, value) in enumerate(initial)
+                    rhs += value * iv^(index - 1) / gamma(index)
+                end
+                for index in range(M, N-1; step=1)
+                    new_z = Symbol(:ʐ, :_, i)
+                    i += 1
+                    γ = γ_i(index)
+                    base = sub_eq
+                    for k in range(1, m; step=1)
+                        new_z = Symbol(:ʐ, :_, index-M, :_, k)
+                        new_z = ModelingToolkit.unwrap(only(@variables $new_z(iv)))
+                        new_eq = D(new_z) ~ base - γ*new_z
+                        base = k * new_z
+                        push!(new_eqs, new_eq)
+                        push!(def, new_z=>0)
+                    end
+                    rhs += coeff*c_i(index)*new_z
+                end
+            end
+        end
+        push!(new_eqs, sub_var ~ rhs)
+        return (new_eqs, def)
+    end
+
+    for (eq, cur_var, alpha, init) in zip(eqs, variables, alphas, initials)
+        (new_eqs, def) = fto_helper(eq, cur_var, alpha; initial=init)
+        append!(all_eqs, new_eqs)
+        append!(all_def, def)
+    end
+    append!(all_eqs, additional_eqs)
+    @named sys = System(all_eqs, iv; defaults=all_def)
+    return mtkcompile(sys)
+end
+
+"""
+Generates the system of ODEs to find solution to FDEs.
+
+Example:
+
+```julia
+@independent_variables t
+@variables x_0(t)
+D = Differential(t)
+tspan = (0., 5000.)
+
+function expect(t)
+    return sqrt(2) * sin(t + pi/4)
+end
+
+sys = linear_fractional_to_ordinary([3, 2.5, 2, 1, .5, 0], [1, 1, 1, 4, 1, 4], 6*cos(t), 10^-5, 5000; initials=[1, 1, -1])
+prob = ODEProblem(sys, [], tspan)
+sol = solve(prob, radau5(), abstol = 1e-5, reltol = 1e-5)
+```
+"""
+function linear_fractional_to_ordinary(degrees, coeffs, rhs, epsilon, T; initials = 0, symbol = :x, iv = only(@independent_variables t))
+    previous = Symbol(symbol, :_, 0)
+    previous = ModelingToolkit.unwrap(only(@variables $previous(iv)))
+    @variables x_0(iv)
+    D = Differential(iv)
+    i = 0
+    all_eqs = Equation[]
+    all_def = Pair{Num, Int64}[]
+
+    function fto_helper(sub_eq, α)
+        δ = (gamma(α+1) * epsilon)^(1/α)
+        a = pi/2*(1-(1-α)/((2-α) * log(epsilon^-1)))
+        h = 2*pi*a / log(1 + (2/epsilon * (cos(a))^(α - 1)))
+
+        x_sub = (gamma(2-α) * epsilon)^(1/(1-α))
+        x_sup = -log(gamma(1-α) * epsilon)
+        M = floor(Int, log(x_sub / T) / h)
+        N = ceil(Int, log(x_sup / δ) / h)
+
+        function c_i(index)
+            h * sin(pi * α) / pi * exp((1-α)*h*index)
+        end
+
+        function γ_i(index)
+            exp(h * index)
+        end
+
+        new_eqs = Equation[]
+        def = Pair{Num, Int64}[]
+        sum = 0
+        for index in range(M, N-1; step=1)
+            new_z = Symbol(:ʐ, :_, i)
+            i += 1
+            new_z = ModelingToolkit.unwrap(only(@variables $new_z(iv)))
+            new_eq = D(new_z) ~ sub_eq - γ_i(index)*new_z
+            push!(new_eqs, new_eq)
+            push!(def, new_z=>0)
+            sum += c_i(index)*new_z
+        end
+        return (new_eqs, def, sum)
+    end
+
+    for i in range(1, ceil(Int, degrees[1]); step=1)
+        new_x = Symbol(symbol, :_, i)
+        new_x = ModelingToolkit.unwrap(only(@variables $new_x(iv)))
+        push!(all_eqs, D(previous) ~ new_x)
+        push!(all_def, previous => initials[i])
+        previous = new_x
+    end
+
+    new_rhs = -rhs
+    for (degree, coeff) in zip(degrees, coeffs)
+        rounded = ceil(Int, degree)
+        new_x = Symbol(symbol, :_, rounded)
+        new_x = ModelingToolkit.unwrap(only(@variables $new_x(iv)))
+        if isinteger(degree)
+            new_rhs += coeff * new_x
+        else
+            (new_eqs, def, sum) = fto_helper(new_x, rounded - degree)
+            append!(all_eqs, new_eqs)
+            append!(all_def, def)
+            new_rhs += coeff * sum
+        end
+    end
+    push!(all_eqs, 0 ~ new_rhs)
+    @named sys = System(all_eqs, iv; defaults=all_def)
+    return mtkcompile(sys)
+end
+
 """
     change_independent_variable(
         sys::System, iv, eqs = [];
diff --git a/test/fractional_to_ordinary.jl b/test/fractional_to_ordinary.jl
new file mode 100644
index 0000000000..2c7dc267f3
--- /dev/null
+++ b/test/fractional_to_ordinary.jl
@@ -0,0 +1,85 @@
+using ModelingToolkit, OrdinaryDiffEq, ODEInterfaceDiffEq, SpecialFunctions
+using Test
+
+# Testing for α < 1
+# Uses example 1 from Section 7 of https://arxiv.org/pdf/2506.04188
+@independent_variables t
+@variables x(t)
+D = Differential(t)
+tspan = (0., 1.)
+timepoint = [0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.]
+
+function expect(t, α)
+    return (3/2*t^(α/2) - t^4)^2
+end
+
+α = 0.5
+eqs = (9*gamma(1 + α)/4) - (3*t^(4 - α/2)*gamma(5 + α/2)/gamma(5 - α/2))
+eqs += (gamma(9)*t^(8 - α)/gamma(9 - α)) + (3/2*t^(α/2)-t^4)^3 - x^(3/2)
+sys = fractional_to_ordinary(eqs, x, α, 10^-7, 1)
+
+prob = ODEProblem(sys, [], tspan)
+sol = solve(prob, radau5(), saveat=timepoint, abstol = 1e-10, reltol = 1e-10)
+
+time = 0
+while(time <= 1)
+    @test isapprox(expect(time, α), sol(time, idxs=x), atol=1e-7)
+    time += 0.1
+end
+
+α = 0.3
+eqs = (9*gamma(1 + α)/4) - (3*t^(4 - α/2)*gamma(5 + α/2)/gamma(5 - α/2))
+eqs += (gamma(9)*t^(8 - α)/gamma(9 - α)) + (3/2*t^(α/2)-t^4)^3 - x^(3/2)
+sys = fractional_to_ordinary(eqs, x, α, 10^-7, 1; matrix=true)
+
+prob = ODEProblem(sys, [], tspan)
+sol = solve(prob, radau5(), saveat=timepoint, abstol = 1e-10, reltol = 1e-10)
+
+time = 0
+while(time <= 1)
+    @test isapprox(expect(time, α), sol(time, idxs=x), atol=1e-7)
+    time += 0.1
+end
+
+α = 0.9
+eqs = (9*gamma(1 + α)/4) - (3*t^(4 - α/2)*gamma(5 + α/2)/gamma(5 - α/2))
+eqs += (gamma(9)*t^(8 - α)/gamma(9 - α)) + (3/2*t^(α/2)-t^4)^3 - x^(3/2)
+sys = fractional_to_ordinary(eqs, x, α, 10^-7, 1)
+
+prob = ODEProblem(sys, [], tspan)
+sol = solve(prob, radau5(), saveat=timepoint, abstol = 1e-10, reltol = 1e-10)
+
+time = 0
+while(time <= 1)
+    @test isapprox(expect(time, α), sol(time, idxs=x), atol=1e-7)
+    time += 0.1
+end
+
+# # Testing for example 2 of Section 7
+# @independent_variables t
+# @variables x(t) y(t)
+# D = Differential(t)
+# tspan = (0., 220.)
+
+# sys = fractional_to_ordinary([1 - 4*x + x^2 * y, 3*x - x^2 * y], [x, y], [1.3, 0.8], 10^-8, 220; initials=[[1.2, 1], 2.8]; matrix=true)
+# prob = ODEProblem(sys, [], tspan)
+# sol = solve(prob, radau5(), abstol = 1e-8, reltol = 1e-8)
+
+# @test isapprox(1.0097684171, sol(220, idxs=x), atol=1e-5)
+# @test isapprox(2.1581264031, sol(220, idxs=y), atol=1e-5)
+
+# #Testing for example 3 of Section 7
+# @independent_variables t
+# @variables x_0(t)
+# D = Differential(t)
+# tspan = (0., 5000.)
+
+# function expect(t)
+#     return sqrt(2) * sin(t + pi/4)
+# end
+
+# sys = linear_fractional_to_ordinary([3, 2.5, 2, 1, .5, 0], [1, 1, 1, 4, 1, 4], 6*cos(t), 10^-5, 5000; initials=[1, 1, -1])
+# prob = ODEProblem(sys, [], tspan)
+# sol = solve(prob, radau5(), abstol = 1e-5, reltol = 1e-5)
+
+# @test isapprox(expect(5000), sol(5000, idxs=x_0), atol=1e-5)
\ No newline at end of file