diff --git a/README.md b/README.md
index d72c4c2b..76397048 100644
--- a/README.md
+++ b/README.md
@@ -107,10 +107,7 @@ By using the ```NonLinearAlg``` algorithms to solve this problem:
 ```julia
 using FractionalDiffEq, Plots
 function chua!(du, x, p, t)
-    a = 10.725; b = 10.593
-    c = 0.268
-    m0 = -1.1726
-    m1 = -0.7872
+    a, b, c, m0, m1 = p
     du[1] = a*(x[2]-x[1]-(m1*x[1]+0.5*(m0-m1)*(abs(x[1]+1)-abs(x[1]-1))))
     du[2] = x[1]-x[2]+x[3]
     du[3] = -b*x[2]-c*x[3]
@@ -118,9 +115,10 @@ end
 α = [0.93, 0.99, 0.92];
 x0 = [0.2; -0.1; 0.1];
 h = 0.01; tspan = (0, 100);
-prob = FODESystem(chua!, α, x0, tspan)
-result = solve(prob, h, NonLinearAlg())
-plot(result[:, 1], result[:, 2], title="Chua System", legend=:bottomright)
+p = [10.725, 10.593, 0.268, -1.1726, -0.7872]
+prob = FODESystem(chua!, α, x0, tspan, p)
+sol = solve(prob, h, NonLinearAlg())
+plot(sol, vars=(1, 2), title="Chua System", legend=:bottomright)
 ```
 
 And plot the result:
diff --git a/docs/src/fdesystem.md b/docs/src/fdesystem.md
index 75275adb..fbf40d14 100644
--- a/docs/src/fdesystem.md
+++ b/docs/src/fdesystem.md
@@ -28,7 +28,7 @@ Use the ```NonLinearAlg``` algorithm in FractionalDiffEq.jl to solve the Chua sy
 ```julia
 using FractionalDiffEq, Plots
 function chua!(du, x, p, t)
-    a = 10.725; b = 10.593; c = 0.268; m0 = -1.1726; m1 = -0.7872;
+    a, b, c, m0, m1 = p
     du[1] = a*(x[2]-x[1]-(m1*x[1]+0.5*(m0-m1)*(abs(x[1]+1)-abs(x[1]-1))))
     du[2] = x[1]-x[2]+x[3]
     du[3] = -b*x[2]-c*x[3]
@@ -36,9 +36,10 @@ end
 α = [0.93, 0.99, 0.92];
 x0 = [0.2; -0.1; 0.1];
 h = 0.01; tspan = (0, 50);
-prob = FODESystem(chua!, α, x0, tspan)
-result = solve(prob, h, NonLinearAlg())
-plot(result[:, 1], result[:, 2], title="Chua System", legend=:bottomright)
+p = [10.725, 10.593, 0.268, -1.1726, -0.7872]
+prob = FODESystem(chua!, α, x0, tspan, p)
+sol = solve(prob, h, NonLinearAlg())
+plot(sol, vars=(1, 2) title="Chua System", legend=:bottomright)
 ```
 
 ![Chua](./assets/chua.png)
diff --git a/examples/chua.jl b/examples/chua.jl
index 228cc5b7..da84ef7b 100644
--- a/examples/chua.jl
+++ b/examples/chua.jl
@@ -1,17 +1,15 @@
 using FractionalDiffEq, Plots
 function chua!(du, x, p, t)
-    a = 10.725; b = 10.593
-    c = 0.268
-    m0 = -1.1726
-    m1 = -0.7872
+    a, b, c, m0, m1 = p
     du[1] = a*(x[2]-x[1]-(m1*x[1]+0.5*(m0-m1)*(abs(x[1]+1)-abs(x[1]-1))))
     du[2] = x[1]-x[2]+x[3]
     du[3] = -b*x[2]-c*x[3]
 end
 α = [0.93, 0.99, 0.92]
 x0 = [0.2; -0.1; 0.1]
-h = 0.001; tspan = (0, 20)
-prob = FODESystem(chua!, α, x0, tspan)
-result = solve(prob, h, NonLinearAlg())
+h = 0.001; tspan = (0, 50)
+p = [10.725, 10.593, 0.268, -1.1726, -0.7872]
+prob = FODESystem(chua!, α, x0, tspan, p)
+sol = solve(prob, h, NonLinearAlg())
 
-plot(result[1, :], result[2, :], title="Chua System", legend=:bottomright)
\ No newline at end of file
+plot(sol, vars=(1, 2), title="Fractional Chua System", legend=:bottomright)
\ No newline at end of file
diff --git a/src/fodesystem/NonLinear.jl b/src/fodesystem/NonLinear.jl
index da2c4c09..384013f6 100644
--- a/src/fodesystem/NonLinear.jl
+++ b/src/fodesystem/NonLinear.jl
@@ -14,6 +14,7 @@ struct NonLinearAlg <: AbstractFDEAlgorithm end
 function solve(prob::FODESystem, h, ::NonLinearAlg, L0=1e10)
     @unpack f, α, u0, tspan, p = prob
     t0 = tspan[1]; T = tspan[2]
+    time = collect(t0:h:T)
     n = length(u0)
     m::Int = round(Int, (T-t0)/h)+1
     g = genfun(1)
@@ -44,7 +45,7 @@ function solve(prob::FODESystem, h, ::NonLinearAlg, L0=1e10)
         z[:, k] = x1 - u0
     end
     result = z + repeat(u0, 1, m)
-    return result
+    return FODESystemSolution(time, result)
 end
 
 """
diff --git a/src/singletermfode/PECE.jl b/src/singletermfode/PECE.jl
index 19b6096b..0f840879 100644
--- a/src/singletermfode/PECE.jl
+++ b/src/singletermfode/PECE.jl
@@ -241,6 +241,10 @@ struct FODESolution <: AbstractFDESolution
     u::AbstractArray
 end
 
+struct FODESystemSolution <: AbstractFDESolution
+    t::AbstractArray
+    u::AbstractArray
+end
 
 struct FDifferenceSolution <: AbstractFDESolution
     t::AbstractArray
diff --git a/src/utils.jl b/src/utils.jl
index 5aa2ec29..3ff98c60 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -113,4 +113,57 @@ Fractional differential equation solutions visulization hooks.
 
 @recipe f(sol::DODESolution) = sol.t, sol.u
 
-@recipe f(sol::FFMODESolution) = sol.t, sol.u
\ No newline at end of file
+@recipe f(sol::FFMODESolution) = sol.t, sol.u
+
+@recipe function f(sol::FODESystemSolution; vars=nothing)
+    if typeof(vars) == Nothing # When vars is not specified, show time versus each variable.
+        l = size(sol.u, 1)
+        for i in 1:l
+            @series begin
+                sol.t, sol.u[i, :]
+            end
+        end
+    else
+        index = Int[]
+        for i in vars
+            append!(index, i)
+        end
+        len = length(index)
+        index0 = findall(x->x==0, index) # Find the 0 index, which is the time index.
+        if length(index0) == 0
+            if len == 3
+                sol.u[index[1], :], sol.u[index[2], :], sol.u[index[3], :]
+            elseif len == 2
+                sol.u[index[1], :], sol.u[index[2], :]
+            else
+                error("Plots variable index is not correct.")
+            end
+        elseif length(index0) == 1
+            newindex = deleteat!(index, index0)
+            newlen = length(newindex)
+            if newlen == 2
+                for i in newindex
+                    @series begin
+                        sol.t, sol.u[i, :]
+                    end
+                end
+            elseif newlen == 1
+                sol.t, sol.u[newindex[1], :]
+            end
+        end
+#=
+        if len == 3 # For 3 dimension
+            sol.u[index[1], :], sol.u[index[2], :], sol.u[index[3], :]
+        elseif len == 2 # For 2 dimension
+            
+            if length(index0) == 0
+                sol.u[index[1], :], sol.u[index[2], :]
+            else
+                sol.t, sol.u[index[index0[1]], :]
+            end
+        else
+            error("Plot variable index is not correct.")
+        end
+=#
+    end
+end
\ No newline at end of file
diff --git a/test/FODETests.jl b/test/FODETests.jl
index 5e5a2674..e4cc168f 100644
--- a/test/FODETests.jl
+++ b/test/FODETests.jl
@@ -134,7 +134,7 @@ end
     prob = FODESystem(chua!, α, x0, tspan)
     result = solve(prob, h, NonLinearAlg())
 
-    @test isapprox(result, [  0.2   0.11749     0.074388  0.0733938  0.117483   0.210073
+    @test isapprox(result.u, [  0.2   0.11749     0.074388  0.0733938  0.117483   0.210073
     -0.1  -0.0590683  -0.018475  0.0192931  0.0534393  0.0844175
      0.1   0.224134    0.282208  0.286636   0.246248   0.168693]; atol=1e-3)
 end