From 74382afc9e0ae63b342501b1e1361428cb7409d3 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Wed, 9 Oct 2024 22:37:10 +0800
Subject: [PATCH] Initial improvement for FODE solvers

Signed-off-by: Qingyu Qu <2283984853@qq.com>
---
 Project.toml                      |   6 +-
 src/FractionalDiffEq.jl           |   6 +-
 src/fode/bdf.jl                   | 175 +++++++++++++++---------------
 src/fode/implicit_pi_rectangle.jl |  62 +++++------
 src/fode/implicit_pi_trapzoid.jl  |  59 ++++------
 src/fode/newton_gregory.jl        | 107 ++++++++----------
 src/fode/pi_pece.jl               |  67 +++++-------
 src/fode/trapezoid.jl             |  92 ++++++----------
 src/types/problem_utils.jl        |  33 ------
 src/utils.jl                      | 168 ----------------------------
 test/fode.jl                      |   1 +
 11 files changed, 253 insertions(+), 523 deletions(-)
 delete mode 100644 src/types/problem_utils.jl

diff --git a/Project.toml b/Project.toml
index 8017b13a..a2c21f6b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,18 +1,20 @@
 name = "FractionalDiffEq"
 uuid = "c7492dd8-6170-483b-af64-cefb6c377d9a"
 authors = ["Qingyu Qu <erikqqy123@gmail.com>"]
-version = "0.4.0"
+version = "0.5.0"
 
 [deps]
 ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
 InvertedIndices = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -23,12 +25,14 @@ ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24"
 [compat]
 ConcreteStructs = "0.2.2"
 DiffEqBase = "6.156"
+FastClosures = "0.3.2"
 FFTW = "1.8.0"
 ForwardDiff = "0.10.36"
 HypergeometricFunctions = "0.3.23"
 InvertedIndices = "1.3"
 Polynomials = "3, 4"
 RecipesBase = "1.3.4"
+RecursiveArrayTools = "3.27"
 Reexport = "1.2.2"
 SciMLBase = "2.55"
 SparseArrays = "1.10"
diff --git a/src/FractionalDiffEq.jl b/src/FractionalDiffEq.jl
index dd22020e..07890c39 100644
--- a/src/FractionalDiffEq.jl
+++ b/src/FractionalDiffEq.jl
@@ -1,7 +1,7 @@
 module FractionalDiffEq
 
-using LinearAlgebra, Reexport, SciMLBase, SpecialFunctions, SparseArrays, ToeplitzMatrices,
-      FFTW, RecipesBase, ForwardDiff, Polynomials, HypergeometricFunctions, DiffEqBase, ConcreteStructs
+using LinearAlgebra, Reexport, SciMLBase, SpecialFunctions, SparseArrays, ToeplitzMatrices, RecursiveArrayTools,
+      FFTW, ForwardDiff, Polynomials, HypergeometricFunctions, DiffEqBase, ConcreteStructs, FastClosures
 
 import SciMLBase: __solve
 import DiffEqBase: solve
@@ -16,8 +16,6 @@ include("types/problems.jl")
 include("types/algorithms.jl")
 include("types/solutions.jl")
 
-include("types/problem_utils.jl")
-
 # Multi-terms fractional ordinary differential equations
 include("multitermsfode/matrix.jl")
 include("multitermsfode/hankelmatrix.jl")
diff --git a/src/fode/bdf.jl b/src/fode/bdf.jl
index ef56e282..27901d09 100644
--- a/src/fode/bdf.jl
+++ b/src/fode/bdf.jl
@@ -8,7 +8,7 @@
     y
     fy
     zn
-    Jfdefun
+    jac
     p
     problem_size
 
@@ -65,12 +65,30 @@ function SciMLBase.__init(prob::FODEProblem, alg::BDF; dt = 0.0, reltol = 1e-6,
     NNr = 2^(Q + 1) * r
 
     # Preallocation of some variables
-    y = zeros(T, problem_size, N + 1)
-    fy = zeros(T, problem_size, N + 1)
+    y = [Vector{T}(undef, problem_size) for _ in 1:(N + 1)]
+    fy = similar(y)
     zn = zeros(T, problem_size, NNr + 1)
 
-    # generate jacobian of the input function
-    Jfdefun(t, u) = jacobian_of_fdefun(prob.f, t, u, p)
+    if prob.f.jac === nothing
+        if iip
+            jac = (df, u, p, t) -> begin
+                _du = similar(u)
+                prob.f(_du, u, p, t)
+                _f = @closure (du, u) -> prob.f(du, u, p, t)
+                ForwardDiff.jacobian!(df, _f, _du, u)
+                return
+            end
+        else
+            jac = (df, u, p, t) -> begin
+                _du = prob.f(u, p, t)
+                _f = @closure (du, u) -> (du .= prob.f(u, p, t))
+                ForwardDiff.jacobian!(df, _f, _du, u)
+                return
+            end
+        end
+    else
+        jac = prob.f.jac
+    end
 
     # Evaluation of convolution and starting BDF_weights of the FLMM
     (omega, w, s) = BDF_weights(alpha, NNr + 1)
@@ -78,19 +96,18 @@ function SciMLBase.__init(prob::FODEProblem, alg::BDF; dt = 0.0, reltol = 1e-6,
 
     # Initializing solution and proces of computation
     mesh = t0 .+ collect(0:N) * dt
-    y[:, 1] = high_order_prob ? u0[1, :] : u0
+    y[1] .= high_order_prob ? u0[1, :] : u0
     temp = high_order_prob ? similar(u0[1, :]) : similar(u0)
-    f(temp, u0, p, t0)
-    fy[:, 1] = temp
+    prob.f(temp, u0, p, t0)
+    fy[1] = temp
 
-    return BDFCache{iip, T}(prob, alg, mesh, u0, alpha, halpha, y, fy, zn, Jfdefun, p,
+    return BDFCache{iip, T}(prob, alg, mesh, u0, alpha, halpha, y, fy, zn, jac, prob.p,
         problem_size, m_alpha, m_alpha_factorial, r, N, Nr, Q, NNr, omega,
         w, s, dt, reltol, abstol, maxiters, high_order_prob, kwargs)
 end
 
 function SciMLBase.solve!(cache::BDFCache{iip, T}) where {iip, T}
     (; prob, alg, mesh, y, r, N, Q, s, dt) = cache
-    tfinal = mesh[end]
 
     BDF_first_approximations(cache)
     BDF_triangolo(cache, s + 1, r - 1, 0)
@@ -98,22 +115,13 @@ function SciMLBase.solve!(cache::BDFCache{iip, T}) where {iip, T}
     # Main process of computation by means of the FFT algorithm
     nx0 = 0
     ny0 = 0
-    ff = zeros(T, 1, 2^(Q + 2), 1)
-    ff[1:2] = [0 2]
+    ff = zeros(T, 1, 2^(Q + 2))
+    ff[1:2] = [0; 2]
     for q in 0:Q
         L = 2^q
         BDF_disegna_blocchi(cache, L, ff, nx0 + L * r, ny0)
         ff[1:(4 * L)] = [ff[1:(2 * L)]; ff[1:(2 * L - 1)]; 4 * L]
     end
-    # Evaluation solution in TFINAL when TFINAL is not in the mesh
-    if tfinal < mesh[N + 1]
-        c = (tfinal - mesh[N]) / dt
-        mesh[N + 1] = tfinal
-        y[:, N + 1] = (1 - c) * y[:, N] + c * y[:, N + 1]
-    end
-    mesh = mesh[1:(N + 1)]
-    y = y[:, 1:(N + 1)]
-    y = collect(Vector{eltype(y)}, eachcol(y))
 
     return DiffEqBase.build_solution(prob, alg, mesh, y)
 end
@@ -122,28 +130,28 @@ function BDF_disegna_blocchi(
         cache::BDFCache{iip, T}, L::P, ff, nx0::P, ny0::P) where {P <: Integer, iip, T}
     (; r, Nr, N) = cache
 
-    nxi::Int = copy(nx0)
-    nxf::Int = copy(nx0 + L * r - 1)
-    nyi::Int = copy(ny0)
-    nyf::Int = copy(ny0 + L * r - 1)
+    nxi::Int = nx0
+    nxf::Int = nx0 + L * r - 1
+    nyi::Int = ny0
+    nyf::Int = ny0 + L * r - 1
     is = 1
-    s_nxi = zeros(T, N)
-    s_nxf = zeros(T, N)
-    s_nyi = zeros(T, N)
-    s_nyf = zeros(T, N)
-    s_nxi[is] = nxi
-    s_nxf[is] = nxf
-    s_nyi[is] = nyi
-    s_nyf[is] = nyf
+    s_nxi = Vector{T}(undef, N)
+    s_nxf = similar(s_nxi)
+    s_nyi = similar(s_nxi)
+    s_nyf = similar(s_nxi)
+    s_nxi[1] = nxi
+    s_nxf[1] = nxf
+    s_nyi[1] = nyi
+    s_nyf[1] = nyf
     i_triangolo = 0
     stop = false
-    while ~stop
+    while !stop
         stop = (nxi + r - 1 == nx0 + L * r - 1) || (nxi + r - 1 >= Nr - 1)
         BDF_quadrato(cache, nxi, nxf, nyi, nyf)
         BDF_triangolo(cache, nxi, nxi + r - 1, nxi)
         i_triangolo = i_triangolo + 1
 
-        if ~stop
+        if !stop
             if nxi + r - 1 == nxf
                 i_Delta = ff[i_triangolo]
                 Delta = i_Delta * r
@@ -178,7 +186,7 @@ function BDF_quadrato(cache::BDFCache, nxi::P, nxf::P, nyi::P, nyf::P) where {P
     funz_beg = nyi + 1
     funz_end = nyf + 1
     vett_coef = omega[(coef_beg + 1):(coef_end + 1)]
-    vett_funz = [cache.fy[:, funz_beg:funz_end] zeros(problem_size, funz_end - funz_beg + 1)]
+    vett_funz = [reduce(hcat, cache.fy[funz_beg:funz_end]) zeros(problem_size, funz_end - funz_beg + 1)]
     zzn = real(fast_conv(vett_coef, vett_funz))
     cache.zn[:, (nxi + 1):(nxf + 1)] = cache.zn[:, (nxi + 1):(nxf + 1)] +
                                        zzn[:, (nxf - nyf):(end - 1)]
@@ -186,73 +194,65 @@ end
 
 function BDF_triangolo(
         cache::BDFCache{iip, T}, nxi::P, nxf::P, j0) where {P <: Integer, iip, T}
-    (; prob, mesh, problem_size, zn, Jfdefun, N, abstol, maxiters, s, w, omega, halpha, p) = cache
+    (; prob, mesh, problem_size, zn, jac, N, abstol, maxiters, s, w, omega, halpha, p) = cache
+    Jfn0 = Matrix{T}(undef, problem_size, problem_size)
     for n in nxi:min(N, nxf)
         n1 = n + 1
         St = ABM_starting_term(cache, mesh[n1])
 
         Phi = zeros(T, problem_size, 1)
         for j in 0:s
-            Phi = Phi + w[j + 1, n1] * cache.fy[:, j + 1]
+            Phi = Phi + w[j + 1, n1] * cache.fy[j + 1]
         end
         for j in j0:(n - 1)
-            Phi = Phi + omega[n - j + 1] * cache.fy[:, j + 1]
+            Phi = Phi + omega[n - j + 1] * cache.fy[j + 1]
         end
         Phi_n = St + halpha * (zn[:, n1] + Phi)
-
-        yn0 = cache.y[:, n]
-        temp = zeros(T, problem_size)
-        prob.f(temp, yn0, p, mesh[n1])
-        fn0 = temp
-        Jfn0 = Jf_vectorfield(mesh[n1], yn0, Jfdefun)
+        yn0 = cache.y[n]
+        fn0 = similar(yn0)
+        prob.f(fn0, yn0, p, mesh[n1])
+        @views jac(Jfn0, yn0, p, mesh[n1])
         Gn0 = yn0 - halpha * omega[1] * fn0 - Phi_n
         stop = false
         it = 0
         yn1 = similar(yn0)
         fn1 = similar(yn0)
-        while ~stop
+        while !stop
             JGn0 = zeros(T, problem_size, problem_size) + I - halpha * omega[1] * Jfn0
-            yn1 = yn0 - JGn0 \ Gn0
+            yn1 = yn0 - vec(JGn0 \ Gn0)
             prob.f(fn1, yn1, p, mesh[n1])
             Gn1 = yn1 - halpha * omega[1] * fn1 - Phi_n
             it = it + 1
 
             stop = (norm(yn1 - yn0, Inf) < abstol) || (norm(Gn1, Inf) < abstol)
-            if it > maxiters && ~stop
+            if it > maxiters && !stop
                 @warn "Non Convergence"
                 stop = true
             end
 
             yn0 = yn1
             Gn0 = Gn1
-            if ~stop
-                Jfn0 = Jf_vectorfield(mesh[n1], yn0, Jfdefun)
+            if !stop
+                @views jac(Jfn0, yn0, p, mesh[n1])
             end
         end
-        cache.y[:, n1] = yn1
-        cache.fy[:, n1] = fn1
+        cache.y[n1] = yn1
+        cache.fy[n1] = fn1
     end
 end
 
 function BDF_first_approximations(cache::BDFCache{iip, T}) where {iip, T}
-    (; prob, mesh, abstol, problem_size, maxiters, s, halpha, omega, w, Jfdefun, p) = cache
+    (; prob, mesh, abstol, problem_size, maxiters, s, halpha, omega, w, jac, p) = cache
 
     Im = zeros(problem_size, problem_size) + I
     Ims = zeros(problem_size * s, problem_size * s) + I
-    Y0 = zeros(T, s * problem_size, 1)
-    F0 = copy(Y0)
-    B0 = copy(Y0)
+    Y0 = VectorOfArray([cache.y[1] for _ in 1:s])
+    F0 = similar(Y0)
+    B0 = similar(Y0)
     for j in 1:s
-        Y0[((j - 1) * problem_size + 1):(j * problem_size), 1] = cache.y[:, 1]
-        temp = zeros(length(cache.y[:, 1]))
-        prob.f(temp, cache.y[:, 1], p, mesh[j + 1])
-        F0[((j - 1) * problem_size + 1):(j * problem_size), 1] = temp
+        prob.f(F0.u[j], cache.y[1], p, mesh[j + 1])
         St = ABM_starting_term(cache, mesh[j + 1])
-        B0[((j - 1) * problem_size + 1):(j * problem_size), 1] = St +
-                                                                 halpha *
-                                                                 (omega[j + 1] +
-                                                                  w[1, j + 1]) *
-                                                                 cache.fy[:, 1]
+        B0.u[j] = St + halpha * (omega[j + 1] + w[1, j + 1]) * cache.fy[1]
     end
     W = zeros(T, s, s)
     for i in 1:s
@@ -265,32 +265,29 @@ function BDF_first_approximations(cache::BDFCache{iip, T}) where {iip, T}
         end
     end
     W = halpha * kron(W, Im)
-    G0 = Y0 - B0 - W * F0
+    G0 = vec(Y0 - B0) - W * vec(F0)
     JF = zeros(T, s * problem_size, s * problem_size)
+    J_temp = Matrix{T}(undef, problem_size, problem_size)
     for j in 1:s
-        JF[((j - 1) * problem_size + 1):(j * problem_size), ((j - 1) * problem_size + 1):(j * problem_size)] = Jf_vectorfield(
-            mesh[j + 1], cache.y[:, 1], Jfdefun)
+        jac(J_temp, cache.y[1], p, mesh[j + 1])
+        JF[((j - 1) * problem_size + 1):(j * problem_size), ((j - 1) * problem_size + 1):(j * problem_size)] .= J_temp
     end
     stop = false
     it = 0
-    F1 = zeros(T, s * problem_size, 1)
+    F1 = similar(F0)
     Y1 = similar(Y0)
-    while ~stop
+    while !stop
         JG = Ims - W * JF
-        Y1 = Y0 - JG \ G0
+        recursive_unflatten!(Y1, vec(Y0) - JG \ G0)
 
         for j in 1:s
-            temp = zeros(T, length(Y1[((j - 1) * problem_size + 1):(j * problem_size), 1]))
-            prob.f(temp, Y1[((j - 1) * problem_size + 1):(j * problem_size), 1],
-                p, mesh[j + 1])
-            F1[((j - 1) * problem_size + 1):(j * problem_size), 1] = temp
+            prob.f(F1.u[j], Y1.u[j], p, mesh[j + 1])
         end
-        G1 = Y1 - B0 - W * F1
+        G1 = vec(Y1 - B0) - W * vec(F1)
 
         it = it + 1
-
         stop = (norm(Y1 - Y0, Inf) < abstol) || (norm(G1, Inf) < abstol)
-        if it > maxiters && ~stop
+        if it > maxiters && !stop
             @warn "Non Convergence"
             stop = 1
         end
@@ -299,15 +296,14 @@ function BDF_first_approximations(cache::BDFCache{iip, T}) where {iip, T}
         G0 = G1
         if ~stop
             for j in 1:s
-                JF[((j - 1) * problem_size + 1):(j * problem_size), ((j - 1) * problem_size + 1):(j * problem_size)] = Jf_vectorfield(
-                    mesh[j + 1],
-                    Y1[((j - 1) * problem_size + 1):(j * problem_size), 1], Jfdefun)
+                jac(J_temp, Y1.u[j], p, mesh[j+1])
+                JF[((j - 1) * problem_size + 1):(j * problem_size), ((j - 1) * problem_size + 1):(j * problem_size)] .= J_temp
             end
         end
     end
     for j in 1:s
-        cache.y[:, j + 1] = Y1[((j - 1) * problem_size + 1):(j * problem_size), 1]
-        cache.fy[:, j + 1] = F1[((j - 1) * problem_size + 1):(j * problem_size), 1]
+        cache.y[j + 1] = Y1.u[j]
+        cache.fy[j + 1] = F1.u[j]
     end
 end
 
@@ -347,7 +343,7 @@ function BDF_weights(alpha, N)
                                     nn .^ (nu + alpha)
         else
             if i == 0
-                nn_nu_alpha[i + 1, :] = zeros(1, N + 1)
+                nn_nu_alpha[i + 1, :] = zeros(N + 1)
             else
                 nn_nu_alpha[i + 1, :] = gamma(nu + 1) / gamma(nu + 1 + alpha) *
                                         nn .^ (nu + alpha)
@@ -369,7 +365,7 @@ function ABM_starting_term(cache::BDFCache{iip, T}, t) where {iip, T}
     (; u0, m_alpha, mesh, m_alpha_factorial, high_order_prob) = cache
     t0 = mesh[1]
     u0 = high_order_prob ? reshape(u0, 1, length(u0)) : u0
-    ys = zeros(size(u0, 1), 1)
+    ys = zeros(size(u0, 1))
     for k in 1:m_alpha
         ys = ys + (t - t0)^(k - 1) / m_alpha_factorial[k] * u0[:, k]
     end
@@ -405,3 +401,12 @@ function _convert_single_term_to_vectorized_prob!(prob::FODEProblem)
         return new_prob
     end
 end
+
+@views function recursive_unflatten!(y::VectorOfArray, x::AbstractArray)
+    i = 0
+    for yᵢ in y
+        copyto!(yᵢ, x[(i + 1):(i + length(yᵢ))])
+        i += length(yᵢ)
+    end
+    return nothing
+end
\ No newline at end of file
diff --git a/src/fode/implicit_pi_rectangle.jl b/src/fode/implicit_pi_rectangle.jl
index 7b422255..d00494a5 100644
--- a/src/fode/implicit_pi_rectangle.jl
+++ b/src/fode/implicit_pi_rectangle.jl
@@ -64,8 +64,8 @@ function SciMLBase.__init(
     NNr = 2^(Qr + 1) * r
 
     # Preallocation of some variables
-    y = zeros(T, problem_size, N + 1)
-    fy = zeros(T, problem_size, N + 1)
+    y = [Vector{T}(undef, problem_size) for _ in 1:(N + 1)]
+    fy = similar(y)
     zn = zeros(problem_size, NNr + 1)
 
     # generate jacobian of the input function
@@ -135,10 +135,10 @@ function SciMLBase.__init(
 
     # Initializing solution and proces of computation
     mesh = t0 .+ collect(0:N) * dt
-    y[:, 1] = high_order_prob ? u0[1, :] : u0
+    y[1] = high_order_prob ? u0[1, :] : u0
     temp = high_order_prob ? similar(u0[1, :]) : similar(u0)
     f(temp, u0, p, t0)
-    fy[:, 1] = temp
+    fy[1] = temp
     return PIRectCache{iip, T}(
         prob, alg, mesh, u0, order, m_alpha, m_alpha_factorial, y, fy,
         p, problem_size, zn, Jfdefun, r, N, Nr, Qr, NNr, bn, halpha1,
@@ -146,7 +146,6 @@ function SciMLBase.__init(
 end
 function SciMLBase.solve!(cache::PIRectCache{iip, T}) where {iip, T}
     (; prob, alg, mesh, u0, y, r, N, Qr) = cache
-    tfinal = mesh[end]
     PIRect_triangolo(cache, 1, r - 1)
 
     # Main process of computation by means of the FFT algorithm
@@ -163,18 +162,7 @@ function SciMLBase.solve!(cache::PIRectCache{iip, T}) where {iip, T}
         ff[card_ff] = 4 * L
     end
 
-    # Evaluation solution in tfinal when tfinal is not in the mesh
-    if tfinal < mesh[N + 1]
-        c = (tfinal - mesh[N]) / dt
-        mesh[N + 1] = tfinal
-        y[:, N + 1] = (1 - c) * y[:, N] + c * y[:, N + 1]
-    end
-
-    mesh = mesh[1:(N + 1)]
-    y = cache.y[:, 1:(N + 1)]
-    u = collect(Vector{eltype(u0)}, eachcol(y))
-
-    return DiffEqBase.build_solution(prob, alg, mesh, u)
+    return DiffEqBase.build_solution(prob, alg, mesh, cache.y)
 end
 
 function PIRect_disegna_blocchi(
@@ -185,15 +173,15 @@ function PIRect_disegna_blocchi(
     nxf::Int = nx0 + L * r - 1
     nyi::Int = ny0
     nyf::Int = ny0 + L * r - 1
-    is = 1
-    s_nxi = zeros(T, N)
-    s_nxf = zeros(T, N)
-    s_nyi = zeros(T, N)
-    s_nyf = zeros(T, N)
-    s_nxi[is] = nxi
-    s_nxf[is] = nxf
-    s_nyi[is] = nyi
-    s_nyf[is] = nyf
+    is::Int = 1
+    s_nxi = Vector{T}(undef, N)
+    s_nxf = similar(s_nxi)
+    s_nyi = similar(s_nxi)
+    s_nyf = similar(s_nxi)
+    s_nxi[1] = nxi
+    s_nxf[1] = nxf
+    s_nyi[1] = nyi
+    s_nyf[1] = nyf
 
     i_triangolo = 0
     stop = false
@@ -242,9 +230,9 @@ function PIRect_quadrato(cache::PIRectCache{iip, T}, nxi::P, nxf::P,
     Nnxf = min(N, nxf)
 
     if nyi == 0
-        vett_funz = [zeros(problem_size, 1) cache.fy[:, (funz_beg + 1):funz_end]]
+        vett_funz = [zeros(problem_size, 1) reduce(hcat, cache.fy[(funz_beg + 1):funz_end])]
     else
-        vett_funz = cache.fy[:, funz_beg:funz_end]
+        vett_funz = reduce(hcat, cache.fy[funz_beg:funz_end])
     end
     vett_funz_fft = rowfft(vett_funz, coef_end)
     zzn = zeros(problem_size, coef_end)
@@ -270,16 +258,16 @@ function PIRect_triangolo(
         n1 = n + 1
         St = PIRect_system_starting_term(cache, mesh[n + 1])
         # Evaluation of the predictor
-        Phi = zeros(problem_size, 1)
+        Phi = zeros(problem_size)
         for j in nxi:(n - 1)
-            Phi = Phi + bn[1:alpha_length, n - j + 1] .* cache.fy[:, j + 1]
+            Phi = Phi + bn[1:alpha_length, n - j + 1] .* cache.fy[j + 1]
         end
         Phi_n = St + halpha1 .* (cache.zn[:, n + 1] + Phi)
         i_alpha_1 = findall(alpha -> abs(alpha - 1) < 1e-14, order)
-        Phi_n[i_alpha_1] = cache.y[i_alpha_1, n]
+        Phi_n[i_alpha_1] = cache.y[n][i_alpha_1]
 
-        yn0 = cache.y[:, n]
-        fn0 = zeros(T, problem_size)
+        yn0 = cache.y[n]
+        fn0 = similar(yn0)
         Jfn0 = zeros(T, problem_size, problem_size)
         prob.f(fn0, yn0, p, mesh[n + 1])
         Jfn0 = Jf_vectorfield(mesh[n + 1], yn0, Jfdefun)
@@ -292,7 +280,7 @@ function PIRect_triangolo(
         while ~stop
             temp = high_order_prob ? diagm([halpha1]) : diagm(halpha1)
             JGn0 = zeros(T, problem_size, problem_size) + I - temp * Jfn0
-            yn1 = yn0 - JGn0 \ Gn0
+            yn1 = yn0 - vec(JGn0 \ Gn0)
             prob.f(fn1, yn1, p, mesh[n + 1])
             Gn1 = yn1 - halpha1 .* fn1 - Phi_n
             it = it + 1
@@ -309,8 +297,8 @@ function PIRect_triangolo(
                 Jfn0 = Jf_vectorfield(mesh[n1], yn0, Jfdefun)
             end
         end
-        cache.y[:, n1] = yn1
-        cache.fy[:, n1] = fn1
+        cache.y[n1] = yn1
+        cache.fy[n1] = fn1
     end
 end
 
@@ -318,7 +306,7 @@ function PIRect_system_starting_term(cache::PIRectCache{iip, T}, t) where {iip,
     (; mesh, u0, m_alpha, m_alpha_factorial, high_order_prob) = cache
     t0 = mesh[1]
     u0 = high_order_prob ? reshape(u0, 1, length(u0)) : u0
-    ys = zeros(size(u0, 1), 1)
+    ys = zeros(size(u0, 1))
     for k in 1:maximum(m_alpha)
         if length(m_alpha) == 1
             ys = ys .+ (t - t0)^(k - 1) / m_alpha_factorial[k] * u0[:, k]
diff --git a/src/fode/implicit_pi_trapzoid.jl b/src/fode/implicit_pi_trapzoid.jl
index 07d4aeef..ff0b355e 100644
--- a/src/fode/implicit_pi_trapzoid.jl
+++ b/src/fode/implicit_pi_trapzoid.jl
@@ -65,8 +65,8 @@ function SciMLBase.__init(
     NNr = 2^(Qr + 1) * r
 
     # Preallocation of some variables
-    y = zeros(T, problem_size, N + 1)
-    fy = zeros(T, problem_size, N + 1)
+    y = [Vector{T}(undef, problem_size) for _ in 1:(N + 1)]
+    fy = similar(y)
     zn = zeros(problem_size, NNr + 1)
 
     # generate jacobian of the input function
@@ -143,10 +143,10 @@ function SciMLBase.__init(
 
     # Initializing solution and proces of computation
     mesh = t0 .+ collect(0:N) * dt
-    y[:, 1] = high_order_prob ? u0[1, :] : u0
+    y[1] = high_order_prob ? u0[1, :] : u0
     temp = high_order_prob ? similar(u0[1, :]) : similar(u0)
     f(temp, u0, p, t0)
-    fy[:, 1] = temp
+    fy[1] = temp
     return PITrapCache{iip, T}(
         prob, alg, mesh, u0, order, m_alpha, m_alpha_factorial, y, fy, p,
         problem_size, zn, Jfdefun, r, N, Nr, Qr, NNr, an, a0, halpha2,
@@ -172,18 +172,7 @@ function SciMLBase.solve!(cache::PITrapCache{iip, T}) where {iip, T}
         ff[card_ff] = 4 * L
     end
 
-    # Evaluation solution in tfinal when tfinal is not in the mesh
-    if tfinal < mesh[N + 1]
-        c = (tfinal - mesh[N]) / dt
-        mesh[N + 1] = tfinal
-        y[:, N + 1] = (1 - c) * y[:, N] + c * y[:, N + 1]
-    end
-
-    mesh = mesh[1:(N + 1)]
-    y = cache.y[:, 1:(N + 1)]
-    u = collect(Vector{eltype(u0)}, eachcol(y))
-
-    return DiffEqBase.build_solution(prob, alg, mesh, u)
+    return DiffEqBase.build_solution(prob, alg, mesh, y)
 end
 
 function PITrap_disegna_blocchi(
@@ -194,15 +183,15 @@ function PITrap_disegna_blocchi(
     nxf::Int = nx0 + L * r - 1
     nyi::Int = ny0
     nyf::Int = ny0 + L * r - 1
-    is = 1
-    s_nxi = zeros(T, N)
-    s_nxf = zeros(T, N)
-    s_nyi = zeros(T, N)
-    s_nyf = zeros(T, N)
-    s_nxi[is] = nxi
-    s_nxf[is] = nxf
-    s_nyi[is] = nyi
-    s_nyf[is] = nyf
+    is::Int = 1
+    s_nxi = Vector{T}(undef, N)
+    s_nxf = similar(s_nxi)
+    s_nyi = similar(s_nxi)
+    s_nyf = similar(s_nxi)
+    s_nxi[1] = nxi
+    s_nxf[1] = nxf
+    s_nyi[1] = nyi
+    s_nyf[1] = nyf
 
     i_triangolo = 0
     stop = false
@@ -251,9 +240,9 @@ function PITrap_quadrato(cache::PITrapCache{iip, T}, nxi::P, nxf::P,
     Nnxf = min(N, nxf)
 
     if nyi == 0
-        vett_funz = [zeros(problem_size, 1) cache.fy[:, (funz_beg + 1):funz_end]]
+        vett_funz = [zeros(problem_size, 1) reduce(hcat, cache.fy[(funz_beg + 1):funz_end])]
     else
-        vett_funz = cache.fy[:, funz_beg:funz_end]
+        vett_funz = reduce(hcat, cache.fy[funz_beg:funz_end])
     end
     vett_funz_fft = rowfft(vett_funz, coef_end)
     zzn = zeros(problem_size, coef_end)
@@ -281,14 +270,14 @@ function PITrap_triangolo(
         # Evaluation of the predictor
         Phi = zeros(problem_size, 1)
         for j in nxi:(n - 1)
-            Phi = Phi + an[1:alpha_length, n - j + 1] .* cache.fy[:, j + 1]
+            Phi = Phi + an[1:alpha_length, n - j + 1] .* cache.fy[j + 1]
         end
         Phi_n = St .+
                 halpha2 .*
-                (a0[1:alpha_length, n + 1] .* cache.fy[:, 1] + cache.zn[:, n + 1] + Phi)
+                (a0[1:alpha_length, n + 1] .* cache.fy[1] + cache.zn[:, n + 1] + Phi)
 
-        yn0 = cache.y[:, n]
-        fn0 = zeros(T, problem_size)
+        yn0 = cache.y[n]
+        fn0 = similar(yn0)
         Jfn0 = zeros(T, problem_size, problem_size)
         prob.f(fn0, yn0, p, mesh[n + 1])
         Jfn0 = Jf_vectorfield(mesh[n + 1], yn0, Jfdefun)
@@ -302,7 +291,7 @@ function PITrap_triangolo(
         while ~stop
             temp = high_order_prob ? diagm([halpha2]) : diagm(halpha2)
             JGn0 = zeros(T, problem_size, problem_size) + I - temp * Jfn0
-            yn1 = yn0 - JGn0 \ Gn0
+            yn1 = yn0 - vec(JGn0 \ Gn0)
             prob.f(fn1, yn1, p, mesh[n + 1])
             Gn1 = yn1 - halpha2 .* fn1 - Phi_n
             it = it + 1
@@ -319,8 +308,8 @@ function PITrap_triangolo(
                 Jfn0 = Jf_vectorfield(mesh[n1], yn0, Jfdefun)
             end
         end
-        cache.y[:, n1] = yn1
-        cache.fy[:, n1] = fn1
+        cache.y[n1] = yn1
+        cache.fy[n1] = fn1
     end
 end
 
@@ -328,7 +317,7 @@ function PITrap_system_starting_term(cache::PITrapCache{iip, T}, t) where {iip,
     (; mesh, u0, m_alpha, m_alpha_factorial, high_order_prob) = cache
     t0 = mesh[1]
     u0 = high_order_prob ? reshape(u0, 1, length(u0)) : u0
-    ys = zeros(size(u0, 1), 1)
+    ys = zeros(size(u0, 1))
     for k in 1:maximum(m_alpha)
         if length(m_alpha) == 1
             ys = ys .+ (t - t0)^(k - 1) / m_alpha_factorial[k] * u0[:, k]
diff --git a/src/fode/newton_gregory.jl b/src/fode/newton_gregory.jl
index 6d737d6f..84418142 100644
--- a/src/fode/newton_gregory.jl
+++ b/src/fode/newton_gregory.jl
@@ -64,8 +64,8 @@ function SciMLBase.__init(prob::FODEProblem, alg::NewtonGregory; dt = 0.0,
     NNr = 2^(Q + 1) * r
 
     # Preallocation of some variables
-    y = zeros(problem_size, N + 1)
-    fy = zeros(problem_size, N + 1)
+    y = [Vector{T}(undef, problem_size) for _ in 1:(N + 1)]
+    fy = similar(y)
     zn = zeros(problem_size, NNr + 1)
 
     # generate jacobian of input function
@@ -77,10 +77,10 @@ function SciMLBase.__init(prob::FODEProblem, alg::NewtonGregory; dt = 0.0,
 
     # Initializing solution and proces of computation
     mesh = t0 .+ collect(0:N) * dt
-    y[:, 1] = high_order_prob ? u0[1, :] : u0
+    y[1] .= high_order_prob ? u0[1, :] : u0
     temp = high_order_prob ? similar(u0[1, :]) : similar(u0)
     f(temp, u0, p, t0)
-    fy[:, 1] = temp
+    fy[1] = temp
 
     return NewtonGregoryCache{iip, T}(
         prob, alg, mesh, u0, alpha, halpha, y, fy, zn, Jfdefun, p,
@@ -88,8 +88,7 @@ function SciMLBase.__init(prob::FODEProblem, alg::NewtonGregory; dt = 0.0,
         w, s, dt, reltol, abstol, maxiters, high_order_prob, kwargs)
 end
 function SciMLBase.solve!(cache::NewtonGregoryCache{iip, T}) where {iip, T}
-    (; prob, alg, mesh, u0, y, r, N, Q, s, dt) = cache
-    tfinal = mesh[end]
+    (; prob, alg, mesh, u0, y, r, N, Q, s) = cache
 
     NG_first_approximations(cache)
     NG_triangolo(cache, s + 1, r - 1, 0)
@@ -104,15 +103,6 @@ function SciMLBase.solve!(cache::NewtonGregoryCache{iip, T}) where {iip, T}
         NG_disegna_blocchi!(cache, L, ff, nx0 + L * r, ny0)
         ff[1:(4 * L)] = [ff[1:(2 * L)]; ff[1:(2 * L - 1)]; 4 * L]
     end
-    # Evaluation solution in TFINAL when TFINAL is not in the mesh
-    if tfinal < mesh[N + 1]
-        c = (tfinal - mesh[N]) / dt
-        mesh[N + 1] = tfinal
-        y[:, N + 1] = (1 - c) * y[:, N] + c * y[:, N + 1]
-    end
-    mesh = mesh[1:(N + 1)]
-    y = y[:, 1:(N + 1)]
-    y = collect(Vector{eltype(u0)}, eachcol(y))
 
     return DiffEqBase.build_solution(prob, alg, mesh, y)
 end
@@ -120,28 +110,28 @@ end
 function NG_disegna_blocchi!(cache::NewtonGregoryCache{iip, T}, L::P, ff,
         nx0::P, ny0) where {P <: Integer, iip, T}
     (; r, Nr, N) = cache
-    nxi::Int = copy(nx0)
-    nxf::Int = copy(nx0 + L * r - 1)
-    nyi::Int = copy(ny0)
-    nyf::Int = copy(ny0 + L * r - 1)
+    nxi::Int = nx0
+    nxf::Int = nx0 + L * r - 1
+    nyi::Int = ny0
+    nyf::Int = ny0 + L * r - 1
     is::Int = 1
-    s_nxi = zeros(T, N)
-    s_nxf = zeros(T, N)
-    s_nyi = zeros(T, N)
-    s_nyf = zeros(T, N)
-    s_nxi[is] = nxi
-    s_nxf[is] = nxf
-    s_nyi[is] = nyi
-    s_nyf[is] = nyf
+    s_nxi = Vector{T}(undef, N)
+    s_nxf = similar(s_nxi)
+    s_nyi = similar(s_nxi)
+    s_nyf = similar(s_nxi)
+    s_nxi[1] = nxi
+    s_nxf[1] = nxf
+    s_nyi[1] = nyi
+    s_nyf[1] = nyf
     i_triangolo = 0
     stop = false
-    while ~stop
+    while !stop
         stop = (nxi + r - 1 == nx0 + L * r - 1) || (nxi + r - 1 >= Nr - 1)
         NG_quadrato(cache, nxi, nxf, nyi, nyf)
         NG_triangolo(cache, nxi, nxi + r - 1, nxi)
         i_triangolo = i_triangolo + 1
 
-        if ~stop
+        if !stop
             if nxi + r - 1 == nxf
                 i_Delta = ff[i_triangolo]
                 Delta = i_Delta * r
@@ -177,7 +167,7 @@ function NG_quadrato(cache::NewtonGregoryCache{iip, T}, nxi::P, nxf::P,
     funz_beg = nyi + 1
     funz_end = nyf + 1
     vett_coef = omega[(coef_beg + 1):(coef_end + 1)]
-    vett_funz = [cache.fy[:, funz_beg:funz_end] zeros(problem_size, funz_end - funz_beg + 1)]
+    vett_funz = [reduce(hcat, cache.fy[funz_beg:funz_end]) zeros(problem_size, funz_end - funz_beg + 1)]
     zzn = real(fast_conv(vett_coef, vett_funz))
     cache.zn[:, (nxi + 1):(nxf + 1)] = cache.zn[:, (nxi + 1):(nxf + 1)] +
                                        zzn[:, (nxf - nyf):(end - 1)]
@@ -192,14 +182,14 @@ function NG_triangolo(
 
         Phi = zeros(problem_size, 1)
         for j in 0:s
-            Phi = Phi + w[j + 1, n1] * cache.fy[:, j + 1]
+            Phi = Phi + w[j + 1, n1] * cache.fy[j + 1]
         end
         for j in j0:(n - 1)
-            Phi = Phi + omega[n - j + 1] * cache.fy[:, j + 1]
+            Phi = Phi + omega[n - j + 1] * cache.fy[j + 1]
         end
         Phi_n = St + halpha * (zn[:, n1] + Phi)
 
-        yn0 = cache.y[:, n]
+        yn0 = cache.y[n]
         temp = zeros(length(yn0))
         prob.f(temp, yn0, p, mesh[n1])
         fn0 = temp
@@ -211,7 +201,7 @@ function NG_triangolo(
         fn1 = similar(yn0)
         while ~stop
             JGn0 = zeros(problem_size, problem_size) + I - halpha * omega[1] * Jfn0
-            yn1 = yn0 - JGn0 \ Gn0
+            yn1 = yn0 - vec(JGn0 \ Gn0)
             prob.f(fn1, yn1, p, mesh[n1])
             Gn1 = yn1 - halpha * omega[1] * fn1 - Phi_n
             it = it + 1
@@ -228,8 +218,8 @@ function NG_triangolo(
                 Jfn0 = Jf_vectorfield(mesh[n1], yn0, Jfdefun)
             end
         end
-        cache.y[:, n1] = yn1
-        cache.fy[:, n1] = fn1
+        cache.y[n1] = yn1
+        cache.fy[n1] = fn1
     end
 end
 
@@ -237,22 +227,15 @@ function NG_first_approximations(cache::NewtonGregoryCache{iip, T}) where {iip,
     (; prob, mesh, abstol, problem_size, maxiters, s, halpha, omega, w, Jfdefun, p) = cache
     Im = zeros(problem_size, problem_size) + I
     Ims = zeros(problem_size * s, problem_size * s) + I
-    Y0 = zeros(s * problem_size, 1)
-    F0 = copy(Y0)
-    B0 = copy(Y0)
+    Y0 = VectorOfArray([cache.y[1] for _ in 1:s])
+    F0 = similar(Y0)
+    B0 = similar(Y0)
     for j in 1:s
-        Y0[((j - 1) * problem_size + 1):(j * problem_size), 1] = cache.y[:, 1]
-        temp = zeros(length(cache.y[:, 1]))
-        prob.f(temp, cache.y[:, 1], p, mesh[j + 1])
-        F0[((j - 1) * problem_size + 1):(j * problem_size), 1] = temp
+        prob.f(F0.u[j], cache.y[1], p, mesh[j + 1])
         St = NG_starting_term(cache, mesh[j + 1])
-        B0[((j - 1) * problem_size + 1):(j * problem_size), 1] = St +
-                                                                 halpha *
-                                                                 (omega[j + 1] +
-                                                                  w[1, j + 1]) *
-                                                                 cache.fy[:, 1]
+        B0.u[j] = St + halpha * (omega[j + 1] + w[1, j + 1]) * cache.fy[1]
     end
-    W = zeros(s, s)
+    W = zeros(T, s, s)
     for i in 1:s
         for j in 1:s
             if i >= j
@@ -263,27 +246,24 @@ function NG_first_approximations(cache::NewtonGregoryCache{iip, T}) where {iip,
         end
     end
     W = halpha * kron(W, Im)
-    G0 = Y0 - B0 - W * F0
-    JF = zeros(s * problem_size, s * problem_size)
+    G0 = vec(Y0 - B0) - W * vec(F0)
+    JF = zeros(T, s * problem_size, s * problem_size)
     for j in 1:s
         JF[((j - 1) * problem_size + 1):(j * problem_size), ((j - 1) * problem_size + 1):(j * problem_size)] = Jf_vectorfield(
-            mesh[j + 1], cache.y[:, 1], Jfdefun)
+            mesh[j + 1], cache.y[1], Jfdefun)
     end
     stop = false
     it = 0
-    F1 = zeros(s * problem_size, 1)
+    F1 = similar(F0)
     Y1 = similar(Y0)
     while ~stop
         JG = Ims - W * JF
-        Y1 = Y0 - JG \ G0
+        recursive_unflatten!(Y1, vec(Y0) - JG \ G0)
 
         for j in 1:s
-            temp = zeros(length(Y1[((j - 1) * problem_size + 1):(j * problem_size), 1]))
-            prob.f(temp, Y1[((j - 1) * problem_size + 1):(j * problem_size), 1],
-                p, mesh[j + 1])
-            F1[((j - 1) * problem_size + 1):(j * problem_size), 1] = temp
+            prob.f(F1.u[j], Y1.u[j], p, mesh[j + 1])
         end
-        G1 = Y1 - B0 - W * F1
+        G1 = vec(Y1 - B0) - W * vec(F1)
 
         it = it + 1
 
@@ -298,14 +278,13 @@ function NG_first_approximations(cache::NewtonGregoryCache{iip, T}) where {iip,
         if ~stop
             for j in 1:s
                 JF[((j - 1) * problem_size + 1):(j * problem_size), ((j - 1) * problem_size + 1):(j * problem_size)] = Jf_vectorfield(
-                    mesh[j + 1],
-                    Y1[((j - 1) * problem_size + 1):(j * problem_size), 1], Jfdefun)
+                    mesh[j + 1], Y1.u[j], Jfdefun)
             end
         end
     end
     for j in 1:s
-        cache.y[:, j + 1] = Y1[((j - 1) * problem_size + 1):(j * problem_size), 1]
-        cache.fy[:, j + 1] = F1[((j - 1) * problem_size + 1):(j * problem_size), 1]
+        cache.y[j + 1] = Y1.u[j]
+        cache.fy[j + 1] = F1.u[j]
     end
 end
 
@@ -363,7 +342,7 @@ function NG_starting_term(cache::NewtonGregoryCache{iip, T}, t) where {iip, T}
     (; u0, m_alpha, mesh, m_alpha_factorial, high_order_prob) = cache
     t0 = mesh[1]
     u0 = high_order_prob ? reshape(u0, 1, length(u0)) : u0
-    ys = zeros(size(u0, 1), 1)
+    ys = zeros(size(u0, 1))
     for k in 1:Int64(m_alpha)
         ys = ys + (t - t0)^(k - 1) / m_alpha_factorial[k] * u0[:, k]
     end
diff --git a/src/fode/pi_pece.jl b/src/fode/pi_pece.jl
index b172487e..ef22a0a9 100644
--- a/src/fode/pi_pece.jl
+++ b/src/fode/pi_pece.jl
@@ -69,8 +69,8 @@ function SciMLBase.__init(prob::FODEProblem, alg::PECE; dt = 0.0, abstol = 1e-6,
     NNr = 2^(Qr + 1) * r
 
     # Preallocation of some variables
-    y = zeros(problem_size, N + 1)
-    fy = zeros(problem_size, N + 1)
+    y = [Vector{T}(undef, problem_size) for _ in 1:(N + 1)]
+    fy = similar(y)
     zn_pred = zeros(problem_size, NNr + 1)
     mu > 0 ? (zn_corr = zeros(problem_size, NNr + 1)) : (zn_corr = 0)
 
@@ -155,10 +155,10 @@ function SciMLBase.__init(prob::FODEProblem, alg::PECE; dt = 0.0, abstol = 1e-6,
 
     # Initializing solution and proces of computation
     mesh = t0 .+ collect(0:N) * dt
-    y[:, 1] = high_order_prob ? u0[1, :] : u0
+    y[1] = high_order_prob ? u0[1, :] : u0
     temp = high_order_prob ? similar(u0[1, :]) : similar(u0)
     f(temp, u0, p, t0)
-    fy[:, 1] = temp
+    fy[1] = temp
     return PECECache{iip, T}(
         prob, alg, mesh, u0, order, m_alpha, m_alpha_factorial, y, fy, p,
         problem_size, zn_pred, zn_corr, r, N, Nr, Qr, NNr, an, bn, a0, halpha1,
@@ -167,7 +167,6 @@ end
 
 function SciMLBase.solve!(cache::PECECache{iip, T}) where {iip, T}
     (; prob, alg, mesh, u0, r, N, Qr) = cache
-    tfinal = mesh[N + 1]
     ABM_triangolo(cache, 1, r - 1)
 
     # Main process of computation by means of the FFT algorithm
@@ -183,19 +182,7 @@ function SciMLBase.solve!(cache::PECECache{iip, T}) where {iip, T}
         card_ff = 2 * card_ff
         ff[card_ff] = 4 * L
     end
-
-    # Evaluation solution in T when T is not in the mesh
-    if tfinal < mesh[N + 1]
-        c = (tfinal - mesh[N]) / dt
-        mesh[N + 1] = tfinal
-        cache.y[:, N + 1] = (1 - c) * cache.y[:, N] + c * cache.y[:, N + 1]
-    end
-
-    mesh = mesh[1:(N + 1)]
-    u = cache.y[:, 1:(N + 1)]
-    u = collect(Vector{eltype(u0)}, eachcol(u))
-
-    return DiffEqBase.build_solution(prob, alg, mesh, u)
+    return DiffEqBase.build_solution(prob, alg, mesh, cache.y)
 end
 
 function DisegnaBlocchi(
@@ -206,14 +193,14 @@ function DisegnaBlocchi(
     nyi::Int = ny0
     nyf::Int = ny0 + L * r - 1
     is = 1
-    s_nxi = zeros(N)
-    s_nxf = zeros(N)
-    s_nyi = zeros(N)
-    s_nyf = zeros(N)
-    s_nxi[is] = nxi
-    s_nxf[is] = nxf
-    s_nyi[is] = nyi
-    s_nyf[is] = nyf
+    s_nxi = Vector{T}(undef, N)
+    s_nxf = similar(s_nxi)
+    s_nyi = similar(s_nxi)
+    s_nyf = similar(s_nxi)
+    s_nxi[1] = nxi
+    s_nxf[1] = nxf
+    s_nyi[1] = nyi
+    s_nyf[1] = nyf
 
     i_triangolo = 0
     stop = false
@@ -264,8 +251,8 @@ function ABM_quadrato(cache::PECECache{iip, T}, nxi::P, nxf::P,
     Nnxf = min(N, nxf)
 
     # Evaluation convolution segment for the predictor
-    vett_funz = cache.fy[:, funz_beg:funz_end]
-    vett_funz_fft = rowfft(vett_funz, coef_end)
+    vett_funz = cache.fy[funz_beg:funz_end]
+    vett_funz_fft = rowfft(reduce(hcat, vett_funz), coef_end)
     zzn_pred = zeros(problem_size, coef_end)
     for i in 1:problem_size
         i_alpha::Int = min(alpha_length, i)
@@ -279,7 +266,7 @@ function ABM_quadrato(cache::PECECache{iip, T}, nxi::P, nxf::P,
     # Evaluation convolution segment for the corrector
     if mu > 0
         if nyi == 0 # Evaluation of the lowest square
-            vett_funz = [zeros(problem_size, 1) cache.fy[:, (funz_beg + 1):funz_end]]
+            vett_funz = [zeros(problem_size, 1) reduce(hcat, cache.fy[(funz_beg + 1):funz_end])]
             vett_funz_fft = rowfft(vett_funz, coef_end)
         end
         zzn_corr = zeros(problem_size, coef_end)
@@ -305,14 +292,14 @@ function ABM_triangolo(
 
     for n in nxi:min(N, nxf)
         # Evaluation of the predictor
-        Phi = zeros(T, problem_size, 1)
+        Phi = zeros(T, problem_size)
         if nxi == 1 # Case of the first triangle
             j_beg = 0
         else # Case of any triangle but not the first
             j_beg = nxi
         end
         for j in j_beg:(n - 1)
-            Phi = Phi .+ bn[1:alpha_length, n - j] .* cache.fy[:, j + 1]
+            Phi = Phi .+ bn[1:alpha_length, n - j] .* cache.fy[j + 1]
         end
         St = starting_term(cache, mesh[n + 1])
         y_pred = St .+ halpha1 .* (zn_pred[:, n + 1] .+ Phi)
@@ -321,17 +308,17 @@ function ABM_triangolo(
 
         # Evaluation of the corrector
         if mu == 0
-            cache.y[:, n + 1] = y_pred
-            cache.fy[:, n + 1] = f_pred
+            cache.y[n + 1] = y_pred
+            cache.fy[n + 1] = f_pred
         else
             j_beg = nxi
-            Phi = zeros(problem_size, 1)
+            Phi = zeros(problem_size)
             for j in j_beg:(n - 1)
-                Phi = Phi + an[1:alpha_length, n - j + 1] .* cache.fy[:, j + 1]
+                Phi = Phi + an[1:alpha_length, n - j + 1] .* cache.fy[j + 1]
             end
             Phi_n = St .+
                     halpha2 .*
-                    (a0[1:alpha_length, n + 1] .* cache.fy[:, 1] .+ zn_corr[:, n + 1] .+
+                    (a0[1:alpha_length, n + 1] .* cache.fy[1] .+ zn_corr[:, n + 1] .+
                      Phi)
             yn0 = y_pred
             fn0 = f_pred
@@ -355,8 +342,8 @@ function ABM_triangolo(
                 yn0 = yn1
                 fn0 = fn1
             end
-            cache.y[:, n + 1] = yn1
-            cache.fy[:, n + 1] = fn1
+            cache.y[n + 1] = yn1
+            cache.fy[n + 1] = fn1
         end
     end
 end
@@ -365,13 +352,13 @@ function starting_term(cache::PECECache{iip, T}, t) where {iip, T}
     (; mesh, m_alpha, u0, m_alpha_factorial, high_order_prob) = cache
     t0 = mesh[1]
     u0 = high_order_prob ? reshape(u0, 1, length(u0)) : u0
-    ys = zeros(size(u0, 1), 1)
+    ys = zeros(size(u0, 1))
     for k in 1:maximum(m_alpha)
         if length(m_alpha) == 1
             ys = ys + (t - t0)^(k - 1) / m_alpha_factorial[k] * u0[:, k]
         else
             i_alpha = findall(x -> x >= k, m_alpha)
-            ys[i_alpha, 1] = ys[i_alpha, 1] +
+            ys[i_alpha] = ys[i_alpha, 1] +
                              (t - t0)^(k - 1) * u0[i_alpha, k] ./
                              m_alpha_factorial[i_alpha, k]
         end
diff --git a/src/fode/trapezoid.jl b/src/fode/trapezoid.jl
index 9a2d5830..9d01105c 100644
--- a/src/fode/trapezoid.jl
+++ b/src/fode/trapezoid.jl
@@ -66,8 +66,8 @@ function SciMLBase.__init(prob::FODEProblem, alg::Trapezoid; dt = 0.0,
     NNr = 2^(Q + 1) * r
 
     # Preallocation of some variables
-    y = zeros(T, problem_size, N + 1)
-    fy = zeros(T, problem_size, N + 1)
+    y = [Vector{T}(undef, problem_size) for _ in 1:(N + 1)]
+    fy = similar(y)
     zn = zeros(T, problem_size, NNr + 1)
 
     # generate jacobian of input function
@@ -79,10 +79,10 @@ function SciMLBase.__init(prob::FODEProblem, alg::Trapezoid; dt = 0.0,
 
     # Initializing solution and proces of computation
     mesh = t0 .+ collect(0:N) * dt
-    y[:, 1] = high_order_prob ? u0[1, :] : u0
+    y[1] = high_order_prob ? u0[1, :] : u0
     temp = high_order_prob ? similar(u0[1, :]) : similar(u0)
     f(temp, u0, p, t0)
-    fy[:, 1] = temp
+    fy[1] = temp
     return TrapezoidCache{iip, T}(prob, alg, mesh, u0, alpha, halpha, y, fy, zn, Jfdefun, p,
         problem_size, m_alpha, m_alpha_factorial, r, N, Nr, Q, NNr, omega,
         w, s, dt, reltol, abstol, maxiters, high_order_prob, kwargs)
@@ -105,15 +105,6 @@ function SciMLBase.solve!(cache::TrapezoidCache{iip, T}) where {iip, T}
         TrapDisegnaBlocchi(cache, L, ff, nx0 + L * r, ny0)
         ff[1:(4 * L)] = [ff[1:(2 * L)]; ff[1:(2 * L - 1)]; 4 * L]
     end
-    # Evaluation solution in TFINAL when TFINAL is not in the mesh
-    if tfinal < mesh[N + 1]
-        c = (tfinal - mesh[N]) / dt
-        mesh[N + 1] = tfinal
-        y[:, N + 1] = (1 - c) * y[:, N] + c * y[:, N + 1]
-    end
-    mesh = mesh[1:(N + 1)]
-    y = y[:, 1:(N + 1)]
-    y = collect(Vector{eltype(u0)}, eachcol(y))
 
     return DiffEqBase.build_solution(prob, alg, mesh, y)
 end
@@ -127,14 +118,14 @@ function TrapDisegnaBlocchi(cache::TrapezoidCache{iip, T}, L::P, ff,
     nyi::Int = copy(ny0)
     nyf::Int = copy(ny0 + L * r - 1)
     is = 1
-    s_nxi = zeros(N)
-    s_nxf = zeros(N)
-    s_nyi = zeros(N)
-    s_nyf = zeros(N)
-    s_nxi[is] = nxi
-    s_nxf[is] = nxf
-    s_nyi[is] = nyi
-    s_nyf[is] = nyf
+    s_nxi = Vector{T}(undef, N)
+    s_nxf = similar(s_nxi)
+    s_nyi = similar(s_nxi)
+    s_nyf = similar(s_nxi)
+    s_nxi[1] = nxi
+    s_nxf[1] = nxf
+    s_nyi[1] = nyi
+    s_nyf[1] = nyf
     i_triangolo = 0
     stop = false
     while ~stop
@@ -179,7 +170,7 @@ function TrapQuadrato(cache::TrapezoidCache{iip, T}, nxi::P, nxf::P,
     funz_beg = nyi + 1
     funz_end = nyf + 1
     vett_coef = omega[(coef_beg + 1):(coef_end + 1)]
-    vett_funz = [cache.fy[:, funz_beg:funz_end] zeros(problem_size, funz_end - funz_beg + 1)]
+    vett_funz = [reduce(hcat, cache.fy[funz_beg:funz_end]) zeros(problem_size, funz_end - funz_beg + 1)]
     zzn = real(fast_conv(vett_coef, vett_funz))
     cache.zn[:, (nxi + 1):(nxf + 1)] = cache.zn[:, (nxi + 1):(nxf + 1)] +
                                        zzn[:, (nxf - nyf):(end - 1)]
@@ -192,16 +183,16 @@ function TrapTriangolo(
         n1::Int = n + 1
         St = TrapStartingTerm(cache, mesh[n1])
 
-        Phi = zeros(problem_size, 1)
+        Phi = zeros(problem_size)
         for j in 0:s
-            Phi = Phi + w[j + 1, n1] * cache.fy[:, j + 1]
+            Phi = Phi + w[j + 1, n1] * cache.fy[j + 1]
         end
         for j in j0:(n - 1)
-            Phi = Phi + omega[n - j + 1] * cache.fy[:, j + 1]
+            Phi = Phi + omega[n - j + 1] * cache.fy[j + 1]
         end
         Phi_n = St + halpha * (zn[:, n1] + Phi)
 
-        yn0 = cache.y[:, n]
+        yn0 = cache.y[n]
         temp = zeros(length(yn0))
         prob.f(temp, yn0, p, mesh[n1])
         fn0 = temp
@@ -213,7 +204,7 @@ function TrapTriangolo(
         fn1 = similar(yn0)
         while ~stop
             JGn0 = zeros(problem_size, problem_size) + I - halpha * omega[1] * Jfn0
-            yn1 = yn0 - JGn0 \ Gn0
+            yn1 = yn0 - vec(JGn0 \ Gn0)
             prob.f(fn1, yn1, p, mesh[n1])
             Gn1 = yn1 - halpha * omega[1] * fn1 - Phi_n
             it = it + 1
@@ -230,8 +221,8 @@ function TrapTriangolo(
                 Jfn0 = Jf_vectorfield(mesh[n1], yn0, Jfdefun)
             end
         end
-        cache.y[:, n1] = yn1
-        cache.fy[:, n1] = fn1
+        cache.y[n1] = yn1
+        cache.fy[n1] = fn1
     end
 end
 
@@ -239,20 +230,13 @@ function TrapFirstApproximations(cache::TrapezoidCache{iip, T}) where {iip, T}
     (; prob, mesh, abstol, problem_size, maxiters, s, halpha, omega, w, Jfdefun, p) = cache
     Im = zeros(problem_size, problem_size) + I
     Ims = zeros(problem_size * s, problem_size * s) + I
-    Y0 = zeros(s * problem_size, 1)
-    F0 = copy(Y0)
-    B0 = copy(Y0)
+    Y0 = VectorOfArray([cache.y[1] for _ in 1:s])
+    F0 = similar(Y0)
+    B0 = similar(Y0)
     for j in 1:s
-        Y0[((j - 1) * problem_size + 1):(j * problem_size), 1] = cache.y[:, 1]
-        temp = zeros(length(cache.y[:, 1]))
-        prob.f(temp, cache.y[:, 1], p, mesh[j + 1])
-        F0[((j - 1) * problem_size + 1):(j * problem_size), 1] = temp
+        prob.f(F0.u[j], cache.y[1], p, mesh[j + 1])
         St = TrapStartingTerm(cache, mesh[j + 1])
-        B0[((j - 1) * problem_size + 1):(j * problem_size), 1] = St +
-                                                                 halpha *
-                                                                 (omega[j + 1] +
-                                                                  w[1, j + 1]) *
-                                                                 cache.fy[:, 1]
+        B0.u[j] = St + halpha * (omega[j + 1] + w[1, j + 1]) * cache.fy[1]
     end
     W = zeros(s, s)
     for i in 1:s
@@ -265,29 +249,26 @@ function TrapFirstApproximations(cache::TrapezoidCache{iip, T}) where {iip, T}
         end
     end
     W = halpha * kron(W, Im)
-    G0 = Y0 - B0 - W * F0
+    G0 = vec(Y0 - B0) - W * vec(F0)
     JF = zeros(s * problem_size, s * problem_size)
     for j in 1:s
         JF[((j - 1) * problem_size + 1):(j * problem_size), ((j - 1) * problem_size + 1):(j * problem_size)] = Jf_vectorfield(
-            mesh[j + 1], cache.y[:, 1], Jfdefun)
+            mesh[j + 1], cache.y[1], Jfdefun)
     end
     stop = false
     it = 0
-    F1 = zeros(s * problem_size, 1)
+    F1 = similar(F0)
+    Y1 = similar(Y0)
     while ~stop
         JG = Ims - W * JF
-        global Y1 = Y0 - JG \ G0
+        recursive_unflatten!(Y1, vec(Y0) - JG \ G0)
 
         for j in 1:s
-            temp = zeros(length(Y1[((j - 1) * problem_size + 1):(j * problem_size), 1]))
-            prob.f(temp, Y1[((j - 1) * problem_size + 1):(j * problem_size), 1],
-                p, mesh[j + 1])
-            F1[((j - 1) * problem_size + 1):(j * problem_size), 1] = temp
+            prob.f(F1.u[j], Y1.u[j], p, mesh[j + 1])
         end
-        G1 = Y1 - B0 - W * F1
+        G1 = vec(Y1 - B0) - W * vec(F1)
 
         it = it + 1
-
         stop = (norm(Y1 - Y0, Inf) < abstol) || (norm(G1, Inf) < abstol)
         if it > maxiters && ~stop
             @warn "Non Convergence"
@@ -299,14 +280,13 @@ function TrapFirstApproximations(cache::TrapezoidCache{iip, T}) where {iip, T}
         if ~stop
             for j in 1:s
                 JF[((j - 1) * problem_size + 1):(j * problem_size), ((j - 1) * problem_size + 1):(j * problem_size)] = Jf_vectorfield(
-                    mesh[j + 1],
-                    Y1[((j - 1) * problem_size + 1):(j * problem_size), 1], Jfdefun)
+                    mesh[j + 1], Y1.u[j], Jfdefun)
             end
         end
     end
     for j in 1:s
-        cache.y[:, j + 1] = Y1[((j - 1) * problem_size + 1):(j * problem_size), 1]
-        cache.fy[:, j + 1] = F1[((j - 1) * problem_size + 1):(j * problem_size), 1]
+        cache.y[j + 1] = Y1.u[j]
+        cache.fy[j + 1] = F1.u[j]
     end
 end
 
@@ -370,7 +350,7 @@ function TrapStartingTerm(cache::TrapezoidCache{iip, T}, t) where {iip, T}
     (; u0, m_alpha, mesh, m_alpha_factorial, high_order_prob) = cache
     t0 = mesh[1]
     u0 = high_order_prob ? reshape(u0, 1, length(u0)) : u0
-    ys = zeros(size(u0, 1), 1)
+    ys = zeros(size(u0, 1))
     for k in 1:m_alpha
         ys = ys + (t - t0)^(k - 1) / m_alpha_factorial[k] * u0[:, k]
     end
diff --git a/src/types/problem_utils.jl b/src/types/problem_utils.jl
deleted file mode 100644
index 1360a45d..00000000
--- a/src/types/problem_utils.jl
+++ /dev/null
@@ -1,33 +0,0 @@
-function Base.summary(io::IO, prob::AbstractFDEProblem)
-    type_color, no_color = SciMLBase.get_colorizers(io)
-    print(io,
-        type_color,
-        nameof(typeof(prob)),
-        no_color,
-        " with uType ",
-        type_color,
-        typeof(prob.u0),
-        no_color,
-        " and tType ",
-        type_color,
-        prob.tspan isa Function ? "Unknown" :
-        (prob.tspan === nothing ? "Nothing" : typeof(prob.tspan[1])),
-        no_color,
-        ". In-place: ",
-        type_color,
-        SciMLBase.isinplace(prob),
-        no_color)
-end
-
-function Base.show(io::IO, mime::MIME"text\plain", A::FDDEProblem)
-    summary(io, A)
-    println(io)
-    print(io, "timespan: ")
-    show(io, mime, A.tspan)
-    println(io)
-    print(io, "order: ")
-    show(io, mime, A.order)
-    println(io)
-    print(io, "u0: ")
-    show(io, mime, A.u0)
-end
diff --git a/src/utils.jl b/src/utils.jl
index bdfc997b..17f556bc 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -65,171 +65,3 @@ function Base.show(io::IO, LE::FOLE)
     printstyled("Timespan:", color = :light_blue)
     printstyled("$(LE.t)")
 end
-
-"""
-Fractional differential equation solutions visulization hooks.
-"""
-@recipe f(sol::FODESolution) = sol.t, sol.u
-
-@recipe f(sol::FDifferenceSolution) = sol.t, sol.u
-
-@recipe f(sol::DODESolution) = sol.t, sol.u
-
-@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
-    end
-end
-
-@recipe function f(sol::FDDESystemSolution; 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
-    end
-end
-
-@recipe function f(sol::FFMODESystemSolution; 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
-    end
-end
-
-@recipe function f(LyapunovExponents::FOLE; vars = nothing)
-    if typeof(vars) == Nothing # When vars is not specified, show time versus each variable.
-        l = size(LyapunovExponents.LE, 1)
-        for i in 1:l
-            @series begin
-                LyapunovExponents.t, LyapunovExponents.LE[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
-                LyapunovExponents.LE[index[1], :], LyapunovExponents.LE[index[2], :],
-                LyapunovExponents.LE[index[3], :]
-            elseif len == 2
-                LyapunovExponents.LE[index[1], :], LyapunovExponents.LE[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
-                        LyapunovExponents.t, LyapunovExponents.LE[i, :]
-                    end
-                end
-            elseif newlen == 1
-                LyapunovExponents.t, LyapunovExponents.LE[newindex[1], :]
-            end
-        end
-    end
-end
diff --git a/test/fode.jl b/test/fode.jl
index fca73b87..561806f0 100644
--- a/test/fode.jl
+++ b/test/fode.jl
@@ -73,6 +73,7 @@ end
     analytical_sol = [[i] for i in temp]
     for alg in [BDF(), Trapezoid(), NewtonGregory(), PITrap(), PECE()]#, PIRect()
         sol = solve(prob, alg, dt = dt)
+        @info "$(alg) method"
         @test norm(analytical_sol .- sol.u) < 1e-4
     end
     for alg in [PIEX()]#, PECE()