diff --git a/Project.toml b/Project.toml
index e4f27642..ec59255a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "TaylorIntegration"
 uuid = "92b13dbe-c966-51a2-8445-caca9f8a7d42"
 repo = "https://github.com/PerezHz/TaylorIntegration.jl.git"
-version = "0.16.0"
+version = "0.16.1"
 
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
diff --git a/src/integrator/taylorinteg.jl b/src/integrator/taylorinteg.jl
index c94905c9..1f7871bb 100644
--- a/src/integrator/taylorinteg.jl
+++ b/src/integrator/taylorinteg.jl
@@ -14,17 +14,11 @@ where the Taylor polynomials associated to the solution will be stored, correspo
 field `:p` in [`TaylorSolution`](@ref). See also [`init_psol`](@ref).
 """
 @inline function set_psol!(::Val{true}, psol::Array{Taylor1{U},1}, nsteps::Int, x::Taylor1{U}) where {U<:Number}
-    @inbounds for k in eachindex(x)
-        TaylorSeries.identity!(psol[nsteps], x, k)
-    end
+    @inbounds psol[nsteps] = deepcopy(x)
     return nothing
 end
 @inline function set_psol!(::Val{true}, psol::Array{Taylor1{U},2}, nsteps::Int, x::Vector{Taylor1{U}}) where {U<:Number}
-    @inbounds for i in eachindex(x)
-        for k in eachindex(x[i])
-            TaylorSeries.identity!(psol[i, nsteps], x[i], k)
-        end
-    end
+    @inbounds psol[:,nsteps] .= deepcopy.(x)
     return nothing
 end
 @inline set_psol!(::Val{false}, args...) = nothing
@@ -43,10 +37,10 @@ associated to the solution will be stored, corresponding to field `:p` in
 [`TaylorSolution`](@ref).
 """
 @inline function init_psol(::Val{true}, xv::Array{U,1}, x::Taylor1{U}) where {U<:Number}
-    return [zero(x) for _ in axes(xv, 1)]
+    return Array{Taylor1{U}}(undef, size(xv, 1)-1)
 end
 @inline function init_psol(::Val{true}, xv::Array{U,2}, x::Array{Taylor1{U},1}) where {U<:Number}
-    return [zero(x[1]) for _ in axes(xv, 1), _ in axes(xv, 2)]
+    return Array{Taylor1{U}}(undef, size(xv, 1), size(xv, 2)-1)
 end
 
 # init_psol