diff --git a/src/lyapunov.jl b/src/lyapunov.jl
index b4cceb1e..d4f377fb 100644
--- a/src/lyapunov.jl
+++ b/src/lyapunov.jl
@@ -14,17 +14,23 @@ Computing fractional order Lyapunov exponent of a fractionl order system.
 }
 ```
 """
-function FOLyapunov(fun, order, t_start, h_norm, t_end, u0, h, out)
+function FOLyapunov(fun, order, t_start, h_norm, t_end, u0, h, out, p)
     if is_all_equal(order)
 
-        LE = commensurate_lyapunov(fun, order, t_start, h_norm, t_end, u0, h, out)
+        LE = commensurate_lyapunov(fun, order, t_start, h_norm, t_end, u0, h, out, p)
     else
-        LE = noncommensurate_lyapunov(fun, order, t_start, h_norm, t_end, u0, h, out)
+        LE = noncommensurate_lyapunov(fun, order, t_start, h_norm, t_end, u0, h, out, p)
     end
     return LE
-
 end
 
+# Dispatch for parameters not specified
+FOLyapunov(fun, order, t_start, h_norm, t_end, u0, h, out) = FOLyapunov(fun, order, t_start, h_norm, t_end, u0, h, out, nothing)
+
+#FOLyapunov(sys::FODESystem, h_norm, h, out) = FOLyapunov(sys.f, sys.α, sys.tspan[1], h_norm, sys.tspan[2], sys.u0, h, out)
+FOLyapunov(sys::FODESystem, h_norm, h, out, p) = FOLyapunov(sys.f, sys.α, sys.tspan[1], h_norm, sys.tspan[2], sys.u0, h, out, p)
+FOLyapunov(sys::FODESystem, h_norm, h, out) = FOLyapunov(sys.f, sys.α, sys.tspan[1], h_norm, sys.tspan[2], sys.u0, h, out, nothing)
+
 @inline function is_all_equal(order::AbstractArray)
     length(order) < 2 && return true
     element1 = order[1]
@@ -34,7 +40,7 @@ end
     return true
 end
 
-function noncommensurate_lyapunov(fun, order, t_start, h_norm, t_end, u0, h, out)# TODO: Generate the Lyapunov exponent plot
+function noncommensurate_lyapunov(fun, order, t_start, h_norm, t_end, u0, h, out, p)# TODO: Generate the Lyapunov exponent plot
     ne::Int = length(u0) # System dimension
 
     tspan = Float64[]
@@ -63,7 +69,7 @@ function noncommensurate_lyapunov(fun, order, t_start, h_norm, t_end, u0, h, out
     t = t_start
     LExp = zeros(ne)
     for it=1:n_it
-        prob = FODESystem(extend_fun, q[:], x[:], (t, t+h_norm))
+        prob = FODESystem(extend_fun, q[:], x[:], (t, t+h_norm), p)
         sol = solve(prob, h, PECE())
         Y = sol.u
         t = t+h_norm
@@ -130,7 +136,7 @@ function jacobian_of_fdefun(f, t, y, p)
     end
 end
 
-function commensurate_lyapunov(fun, order, t_start, h_norm, t_end, u0, h, out)# TODO: Generate the Lyapunov exponent plot
+function commensurate_lyapunov(fun, order, t_start, h_norm, t_end, u0, h, out, p)# TODO: Generate the Lyapunov exponent plot
     ne::Int = length(u0) # System dimension
     order = order[1] # Since this is the commensurate case, we only need one element in order array
     
@@ -221,9 +227,6 @@ function commensurate_lyapunov(fun, order, t_start, h_norm, t_end, u0, h, out)#
     return FOLE(tspan, LE)
 end
 
-#FOLyapunov(sys::FODESystem, h_norm, h, out) = FOLyapunov(sys.f, sys.α, sys.tspan[1], h_norm, sys.tspan[2], sys.u0, h, out)
-FOLyapunov(sys::FODESystem, h_norm, h, out) = FOLyapunov(sys.f, sys.α, sys.tspan[1], h_norm, sys.tspan[2], sys.u0, h, out)
-
 mutable struct M
     an
     bn
diff --git a/test/FOLEtest.jl b/test/FOLEtest.jl
index 1f4f307b..608bbac4 100644
--- a/test/FOLEtest.jl
+++ b/test/FOLEtest.jl
@@ -27,4 +27,20 @@ end
     @test isapprox(LE.LE[end-2:end], [-0.0033493403496077058
     -0.0050780674255360226
     -1.774132087623718]; atol=1e-3)
+end
+
+@testset "Test Noncommensurate Lyapunov exponents with parameters specification" begin
+    function LE_RF_TEST(du, u, p, t)
+        a, b, c = p
+        du[1] = u[2]*(u[3]-a+u[1]^2) + b*u[1]
+        du[2] = u[1]*(3*u[3]+1-u[1]^2) + 0.1*u[2]
+        du[3] = -2*u[3]*(c+u[1]*u[2])
+    end
+
+    p = [1, 0.1, 0.98]
+
+    LE = FOLyapunov(LE_RF_TEST, [0.995, 0.992, 0.996], 0, 0.1, 1000, [1,1,1], 0.01, 1000, p)   
+    @test isapprox(LE.LE[end-2:end], [-0.0033493403496077058
+    -0.0050780674255360226
+    -1.774132087623718]; atol=1e-3)
 end
\ No newline at end of file