diff --git a/Project.toml b/Project.toml
index 07611ad0..64a49478 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,6 +6,7 @@ version = "0.3.3"
 [deps]
 ApproxFun = "28f2ccd6-bb30-5033-b560-165f7b14dc2f"
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
+ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
@@ -27,6 +28,7 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 [compat]
 ApproxFun = "0.13"
 ArrayInterface = "7"
+ConcreteStructs = "0.2.2"
 DiffEqBase = "6"
 FFTW = "1.0.0"
 ForwardDiff = "0.10"
diff --git a/README.md b/README.md
index 0676fc00..ed3a98e1 100644
--- a/README.md
+++ b/README.md
@@ -62,9 +62,9 @@ So we can use FractionalDiffEq.jl to solve the problem:
 ```julia
 using FractionalDiffEq, Plots
 fun(u, p, t) = 1-u
-u0 = [0, 0]; tspan = (0, 20); h = 0.001;
+u0 = [0, 0]; tspan = (0, 20);
 prob = SingleTermFODEProblem(fun, 1.8, u0, tspan)
-sol = solve(prob, h, PECE())
+sol = solve(prob, PECE(), dt=0.001)
 plot(sol)
 ```
 
@@ -74,7 +74,7 @@ And if you plot the result, you can see the result of the above IVP:
 
 ### A sophisticated example
 
-Let's see if the initial value problem like:
+Let's see if we have an initial value problem:
 
 $$ y'''(t)+\frac{1}{16}{^C_0D^{2.5}_t}y(t)+\frac{4}{5}y''(t)+\frac{3}{2}y'(t)+\frac{1}{25}{^C_0D^{0.5}_t}y(t)+\frac{6}{5}y(t)=\frac{172}{125}\cos(\frac{4t}{5}) $$
 
@@ -82,10 +82,10 @@ $$ y(0)=0,\ y'(0)=0,\ y''(0)=0 $$
 
 ```julia
 using FractionalDiffEq, Plots
-h=0.01; tspan = (0, 30)
+tspan = (0, 30)
 rightfun(x, y) = 172/125*cos(4/5*x)
 prob = MultiTermsFODEProblem([1, 1/16, 4/5, 3/2, 1/25, 6/5], [3, 2.5, 2, 1, 0.5, 0], rightfun, [0, 0, 0, 0, 0, 0], tspan)
-sol = solve(prob, h, PIEX())
+sol = solve(prob, PIEX(), dt=0.01)
 plot(sol, legend=:bottomright)
 ```
 
@@ -115,10 +115,10 @@ function chua!(du, x, p, t)
 end
 α = [0.93, 0.99, 0.92];
 x0 = [0.2; -0.1; 0.1];
-h = 0.01; tspan = (0, 100);
+tspan = (0, 100);
 p = [10.725, 10.593, 0.268, -1.1726, -0.7872]
 prob = FODESystem(chua!, α, x0, tspan, p)
-sol = solve(prob, h, NonLinearAlg())
+sol = solve(prob, NonLinearAlg(), dt=0.01)
 plot(sol, vars=(1, 2), title="Chua System", legend=:bottomright)
 ```
 
@@ -143,33 +143,14 @@ $$ y(t)=19,\ t<0 $$
 using FractionalDiffEq, Plots
 ϕ(x) = x == 0 ? (return 19.00001) : (return 19.0)
 f(t, y, ϕ) = 3.5*y*(1-ϕ/19)
-h = 0.05; α = 0.97; τ = 0.8; T = 56
-fddeprob = FDDEProblem(f, ϕ, α, τ, T)
-V, y = solve(fddeprob, h, DelayPECE())
-plot(y, V, xlabel="y(t)", ylabel="y(t-τ)")
+α = 0.97; τ = 0.8; tspan = (0,56)
+fddeprob = FDDEProblem(f, ϕ, α, constant_lags=[τ], tspan)
+sol = solve(fddeprob, DelayPECE(), dt=0.05)
+plot(sol, xlabel="y(t)", ylabel="y(t-τ)")
 ```
 
 ![Delayed](docs/src/assets/fdde_example.png)
 
-## Lyapunov exponents of fractional order system
-
-FractionalDiffEq.jl is capable of generating lyapunov exponents of a fractional order system:
-
-Rabinovich-Fabrikant system:
-
-$$
-\begin{cases} D^{\alpha_1} x=y(z-1+z^2)+\gamma x\\
-D^{\alpha_2} y=x(3z+1-x^2)+\gamma y\\
-D^{\alpha_3} z=-2z(\alpha+xy)
-\end{cases}
-$$
-
-```julia
-julia>LE, tspan = FOLyapunov(RF, 0.98, 0, 0.02, 300, [0.1; 0.1; 0.1], 0.005, 1000)
-```
-
-![RF](docs/src/assets/RFLE.png)
-
 # Available Solvers
 
 For more performant solvers, please refer to the [FractionalDiffEq.jl Solvers](https://scifracx.org/FractionalDiffEq.jl/dev/algorithms/) page.
diff --git a/src/FractionalDiffEq.jl b/src/FractionalDiffEq.jl
index 16d6cced..ea8d1e41 100644
--- a/src/FractionalDiffEq.jl
+++ b/src/FractionalDiffEq.jl
@@ -6,6 +6,7 @@ using LinearAlgebra, Reexport, SciMLBase, SpecialFunctions, SparseArrays, Toepli
 
 import DiffEqBase: solve
 import SciMLBase: __solve
+import ConcreteStructs: @concrete
 import InvertedIndices: Not
 import SpecialMatrices: Vandermonde
 import FFTW: fft, ifft
diff --git a/src/fodesystem/bdf.jl b/src/fodesystem/bdf.jl
index e691a00b..13e27d53 100644
--- a/src/fodesystem/bdf.jl
+++ b/src/fodesystem/bdf.jl
@@ -1,3 +1,12 @@
+#=
+@concrete struct BDFCache
+    prob
+    alg
+    mesh
+    
+end
+=#
+
 function SciMLBase.__solve(prob::FODEProblem, alg::BDF; dt=0.0, reltol=1e-6, abstol=1e-6, maxiters = 100)
     dt ≤ 0 ? throw(ArgumentError("dt must be positive")) : nothing
     @unpack f, order, u0, tspan, p = prob