diff --git a/Project.toml b/Project.toml
index e7d780bb2..35fc41cd4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -26,13 +26,13 @@ BenchmarkTools = "1.3"
 DiffEqBase = "6"
 Elliptic = "0.5, 1.0"
 Espresso = "0.6.1"
-OrdinaryDiffEq = "5, 6"
+OrdinaryDiffEq = "6"
 PkgBenchmark = "0.2"
 RecursiveArrayTools = "2"
 Reexport = "0.2, 1"
 Requires = "0.5.2, 1"
 StaticArrays = "0.12.5, 1"
-TaylorSeries = "0.12"
+TaylorSeries = "0.13"
 julia = "1"
 
 [extras]
diff --git a/benchmark/kepler_benchmarks.jl b/benchmark/kepler_benchmarks.jl
index 379d77fce..0d3575290 100644
--- a/benchmark/kepler_benchmarks.jl
+++ b/benchmark/kepler_benchmarks.jl
@@ -115,23 +115,23 @@ let
     SUITE["Kepler"]["kepler6-3"] = @benchmarkable taylorinteg(
         kepler6!, $q0, $t0, $tf3, $_order, $_abstol, -1.0, maxsteps=$maxsteps3)
 
-    # # ==========
-    # # KeplerNotParsed
-    # # ==========
-    # SUITE["KeplerNotParsed"] = BenchmarkGroup()
-
-    # SUITE["KeplerNotParsed"]["kepler1"] = @benchmarkable taylorinteg(
-    #     kepler1!, $q0, $t0, $tf, $_order, $_abstol, -1.0, maxsteps=$maxsteps,
-    #     parse_eqs=false)
-    # SUITE["KeplerNotParsed"]["kepler2"] = @benchmarkable taylorinteg(
-    #     kepler2!, $q0, $t0, $tf, $_order, $_abstol, $pars, maxsteps=$maxsteps,
-    #     parse_eqs=false)
-    #
-    # SUITE["KeplerNotParsed"]["kepler5"] = @benchmarkable taylorinteg(
-    #     kepler5!, $q0, $t0, $tf, $_order, $_abstol, maxsteps=$maxsteps,
-    #     parse_eqs=false)
-    # SUITE["KeplerNotParsed"]["kepler6"] = @benchmarkable taylorinteg(
-    #     kepler6!, $q0, $t0, $tf, $_order, $_abstol, -1.0, maxsteps=$maxsteps,
-    #     parse_eqs=false)
-    #
+    # ==========
+    # KeplerNotParsed
+    # ==========
+    SUITE["KeplerNotParsed"] = BenchmarkGroup()
+
+    SUITE["KeplerNotParsed"]["kepler1"] = @benchmarkable taylorinteg(
+        kepler1!, $q0, $t0, $tf2, $_order, $_abstol, -1.0, maxsteps=$maxsteps2,
+        parse_eqs=false)
+    SUITE["KeplerNotParsed"]["kepler2"] = @benchmarkable taylorinteg(
+        kepler2!, $q0, $t0, $tf2, $_order, $_abstol, $pars, maxsteps=$maxsteps2,
+        parse_eqs=false)
+
+    SUITE["KeplerNotParsed"]["kepler5"] = @benchmarkable taylorinteg(
+        kepler5!, $q0, $t0, $tf2, $_order, $_abstol, maxsteps=$maxsteps2,
+        parse_eqs=false)
+    SUITE["KeplerNotParsed"]["kepler6"] = @benchmarkable taylorinteg(
+        kepler6!, $q0, $t0, $tf2, $_order, $_abstol, -1.0, maxsteps=$maxsteps2,
+        parse_eqs=false)
+
 end
diff --git a/docs/src/taylorize.md b/docs/src/taylorize.md
index 962e3e023..dd43fad90 100644
--- a/docs/src/taylorize.md
+++ b/docs/src/taylorize.md
@@ -17,19 +17,22 @@ order to compute the recurrence relations that are used to construct
 the Taylor expansion of the solution. This is done for each order of
 the series in [`TaylorIntegration.jetcoeffs!`](@ref). These computations are
 not optimized: they waste memory due to repeated allocations of some temporary
-arrays, and perform some operations whose result has been previously
+arrays, and perform some operations whose result has already been previously
 computed.
 
 Here we describe one way to optimize this: The idea is to replace the
 default method of [`TaylorIntegration.jetcoeffs!`](@ref) by another
 method (same function name) which is called by dispatch, and that in principle
 performs better.
-The new method is constructed specifically for the specific function
+The new method is constructed specifically for the function
 defining the equations of motion by parsing its expression. This new
 method performs in principle *exactly* the same operations, but avoids
-the extra allocations; to do the latter, the macro also creates an *internal*
+repeating some operations and the extra allocations.
+To achieve the latter, the macro also creates an *internal*
 function `TaylorIntegration._allocate_jetcoeffs!`, which allocates all temporary
-`Taylor1` objects as well as the declared `Array{Taylor1,1}`s.
+`Taylor1` objects as well as the declared `Array{Taylor1,N}`s, which are stored
+in a [`RetAlloc{T}`](@ref) struct for efficiency, and include arrays (of `Taylor1{T}`
+objects) with up-to-three indices.
 
 
 ## An example
@@ -41,7 +44,7 @@ using the default method, as described [before](@ref pendulum).
 ```@example taylorize
 using TaylorIntegration
 
-function pendulum!(dx, x, p, t)
+function pendulumNP!(dx, x, p, t) # `pendulum!` ODEs, not parsed
     dx[1] = x[2]
     dx[2] = -sin(x[1])
     return dx
@@ -53,28 +56,30 @@ tf = 10000.0
 q0 = [1.3, 0.0]
 
 # The actual integration
-t1, x1 = taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=1500); # warm-up run
-e1 = @elapsed taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=1500);
-all1 = @allocated taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=1500);
+t1, x1 = taylorinteg(pendulumNP!, q0, t0, tf, 25, 1e-20, maxsteps=50000); # warm-up run
+e1 = @elapsed taylorinteg(pendulumNP!, q0, t0, tf, 25, 1e-20, maxsteps=50000);
+all1 = @allocated taylorinteg(pendulumNP!, q0, t0, tf, 25, 1e-20, maxsteps=50000);
 e1, all1
 ```
 
-We note that the initial number of methods defined for
-`TaylorIntegration.jetcoeffs!` is 2.
+The initial number of methods defined for
+`TaylorIntegration.jetcoeffs!` is 2; yet, since `@taylorize` was used
+in [an example previously](@ref diffeqinterface), the current number
+of methods is 3, as explained below.
 ```@example taylorize
-length(methods(TaylorIntegration.jetcoeffs!)) == 2 # initial value
+println(methods(TaylorIntegration.jetcoeffs!)) # default methods
 ```
-Similarly, the number of methods for `TaylorIntegration._allocate_jetcoeffs!` is
-also 2.
+Similarly, the number of methods for `TaylorIntegration._allocate_jetcoeffs!` originally is 2, and for the same reasons it is
+currently 3.
 ```@example taylorize
-length(methods(TaylorIntegration._allocate_jetcoeffs!)) == 2 # initial value
+println(methods(TaylorIntegration._allocate_jetcoeffs!)) # default methods
 ```
-Using `@taylorize` will increase this number by creating a new method for these
-two functions.
+Using `@taylorize` will increase this number by creating a new method
+for these functions.
 
 The macro [`@taylorize`](@ref) is intended to be used in front of the function
 that implements the equations of motion. The macro does the following: it
-first parses the function as it is, so the integration can be computed
+first parses the function as it is, so the integration can still be computed
 using [`taylorinteg`](@ref) as above, by explicitly using the keyword
 argument `parse_eqs=false`; this also declares the function of the ODEs, whose name
 is used for parsing. It then creates and evaluates a new method of
@@ -93,30 +98,30 @@ println(methods(pendulum!))
 ```
 
 ```@example taylorize
-println(methods(TaylorIntegration.jetcoeffs!))
+println(methods(TaylorIntegration.jetcoeffs!)) # result should be 4
 ```
 
 We see that there is only one method of `pendulum!`, and
-there is a *new* method (three in total) of `TaylorIntegration.jetcoeffs!`,
+there is a *new* method (four in total) of `TaylorIntegration.jetcoeffs!`,
 whose signature appears
-in this documentation as `Val{Main.ex-taylorize.pendulum!}`. It is
+in this documentation as `Val{Main.pendulum!}`. It is
 a specialized version for the function `pendulum!` (with some extra information
 about the module where the function was created). This method
 is selected internally if it exists (default), exploiting dispatch, when
-calling [`taylorinteg`](@ref) or [`lyap_taylorinteg`](@ref); to integrate
-using the hard-coded (default) method of [`TaylorIntegration.jetcoeffs!`](@ref) of the
+calling [`taylorinteg`](@ref) or [`lyap_taylorinteg`](@ref). In order to integrate
+using the hard-coded standard (default) method of [`TaylorIntegration.jetcoeffs!`](@ref) of the
 integration above, the keyword argument `parse_eqs` has to be set to `false`.
 Similarly, one can check that there exists a new method of
-`TaylorIntegration._allocate_jetcoeffs!`
+`TaylorIntegration._allocate_jetcoeffs!`.
 
 Now we carry out the integration using the specialized method; note that we
 use the same instruction as above; the default value for the keyword argument `parse_eqs`
 is `true`, so we may omit it.
 
 ```@example taylorize
-t2, x2 = taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=1500); # warm-up run
-e2 = @elapsed taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=1500);
-all2 = @allocated taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=1500);
+t2, x2 = taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=50000); # warm-up run
+e2 = @elapsed taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=50000);
+all2 = @allocated taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=50000);
 e2, all2
 ```
 
@@ -135,10 +140,10 @@ created by [`@taylorize`](@ref), [`taylorinteg`](@ref) and
 [`lyap_taylorinteg`](@ref) recognize the keyword argument `parse_eqs`;
 setting it to `false` causes the standard method to be used.
 ```@example taylorize
-taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=1500, parse_eqs=false); # warm-up run
+taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=50000, parse_eqs=false); # warm-up run
 
-e3 = @elapsed taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=1500, parse_eqs=false);
-all3 = @allocated taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=1500, parse_eqs=false);
+e3 = @elapsed taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=50000, parse_eqs=false);
+all3 = @allocated taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=50000, parse_eqs=false);
 e1/e3, all1/all3
 ```
 
@@ -153,20 +158,21 @@ e4 = @elapsed solve(prob, TaylorMethod(25), abstol=1e-20, parse_eqs=true);
 
 e1/e4
 ```
-Note that there is an additional marginal cost to using `solve` in comparison
-with `taylorinteg`.
+Note that there is an additional cost to using `solve` in comparison
+with `taylorinteg`, but still `@taylorize` yields improved running times.
 
 The speed-up obtained comes from the design of the new (specialized) method of
 `TaylorIntegration.jetcoeffs!` as described [above](@ref idea): it avoids some
 allocations and some repeated computations. This is achieved by knowing the
 specific AST of the function of the ODEs integrated, which is walked
 through and *translated* into the actual implementation, where some
-required auxiliary arrays are created, and the low-level functions defined in
+required auxiliary arrays are created and reused, and the low-level functions defined in
 `TaylorSeries.jl` are used.
 For this, we heavily rely on [`Espresso.jl`](https://github.com/dfdx/Espresso.jl) and
 some metaprogramming; we thank Andrei Zhabinski for his help and comments.
 
-The new `TaylorIntegration.jetcoeffs!` and `TaylorIntegration._allocate_jetcoeffs!` method can be inspected by
+The new `TaylorIntegration.jetcoeffs!` and `TaylorIntegration._allocate_jetcoeffs!`
+methods can be inspected by
 constructing the expression corresponding to the function, and using
 [`TaylorIntegration._make_parsed_jetcoeffs`](@ref):
 
@@ -182,7 +188,7 @@ new_ex1, new_ex2 = TaylorIntegration._make_parsed_jetcoeffs(ex)
 
 The first function has a similar structure as the hard-coded method of
 `TaylorIntegration.jetcoeffs!`, but uses low-level functions in `TaylorSeries`
-(e.g., `sincos!` above); all temporary `Taylor1` objects as well as declared
+(e.g., `sincos!` above). Temporary `Taylor1` objects as well as declared
 arrays are allocated once by `TaylorIntegration._allocate_jetcoeffs!`.
 More complex functions quickly become very difficult to read. Note that,
 if necessary, one can further optimize `new_ex` manually.
@@ -205,13 +211,17 @@ list some limitations and provide some advice.
   example, the expression `x += y` is not recognized by `@taylorize`. Likewise,
   expressions such as `x = x+y` are not supported by `@taylorize` and should be
   replaced by equivalent expressions in which variables appear only on one side
-  of the assignment; e.g. `z = x+y; x = z`.
+  of the assignment; e.g. `z = x+y; x = z`. The introduction of such temporary
+  variables `z` is left to the user.
 
-- The macro allows the use of array declarations through `Array`, but other ways
-  (e.g. `similar`) are not yet implemented.
+- The macro allows the use of array declarations through `Array` or `Vector`, but other ways
+  (e.g. `similar`) are not yet implemented. Note that certain temporary arrays
+  may be introduced to avoid re-computating certain expressions; only up-to-three
+  indices expressions are currently handled.
 
-- Avoid using variables prefixed by an underscore, in particular `_T`, `_S` and
-  `_N`; using them may lead to name collisions with some internal variables.
+- Avoid using variables prefixed by an underscore, in particular `_T`, `_S`,
+  `_N` and `__idx`, as well as `ord`; using them may lead to name collisions
+  with some internal variables used in the constructed expressions.
 
 - Broadcasting is not recognized by `@taylorize`.
 
@@ -226,25 +236,26 @@ list some limitations and provide some advice.
   rather than comparing against numeric literals.
 
 - Input and output lengths should be determined at the time of `@taylorize`
-  application, not at runtime.  Do not use the length of the input as an
+  application (parse time), not at runtime.  Avoid using the length of the input as an
   implicit indicator of whether to write all elements of the output.  If
-  conditional output of auxiliary equations is desired, use explicit methods,
-  such as through parameters or by setting auxiliary `t0` vector elements
+  conditional output of auxiliary equations is desired use explicit methods,
+  such as through parameters or by setting auxiliary vector elements
   to zero, and assigning unwanted auxiliary outputs zero.
 
 - Expressions which correspond to function calls (so the `head` field is
   `:call`) which are not recognized by the parser are simply copied. The
   heuristics used, especially for vectors, may not work for all cases.
 
-- Use `local` for internal parameters (simple constant values); this improves
-  performance. Do not use it if the variable is Taylor expanded during the integration.
-  
+- Use `local` for internal parameters, e.g., simple constant values; this improves
+  performance. Do not use it if the variable is needed to be Taylor expanded
+  during the integration step.
+
 - To examine the code generated for `jetcoeffs!` and `_allocate_jetcoeffs!`
   for a specific ODE function, follow the pendulum example above; create an expression
-  by wrapping the ODE function (without `@taylorize` prefix) in `:()` and
+  by wrapping the ODE function (without `@taylorize` prefix) in a `:()`-block, and
   supply the expression to `TaylorIntegration._make_parsed_jetcoeffs`.  This
   can help in debugging issues with either function generated by `@taylorize`.
-  
+
 - `@taylorize` supports multi-threading via `Threads.@threads`. **WARNING**:
   this feature is experimental. Since thread-safety depends on the definition
   of each ODE, we cannot guarantee the resulting code to be thread-safe in
@@ -252,7 +263,8 @@ list some limitations and provide some advice.
   thread-safe. For more information about multi-threading, the reader is
   referred to the [Julia documentation](https://docs.julialang.org/en/v1/manual/multi-threading/#man-multithreading).
 
-It is recommended to skim `test/taylorize.jl`, which implements different
-cases and highlights cases where the macro doesn't work and how to solve the problem.
+We recommend to skim `test/taylorize.jl`, which implements different
+cases and highlights examples where the macro does not work, and how to solve the problem;
+read the information that is in the comments.
 
 Please report any problems you may encounter.
diff --git a/src/TaylorIntegration.jl b/src/TaylorIntegration.jl
index 334cd6ac8..6fa9a7722 100644
--- a/src/TaylorIntegration.jl
+++ b/src/TaylorIntegration.jl
@@ -12,10 +12,10 @@ using InteractiveUtils: methodswith
 
 export taylorinteg, lyap_taylorinteg, @taylorize
 
+include("parse_eqs.jl")
 include("explicitode.jl")
 include("lyapunovspectrum.jl")
 include("rootfinding.jl")
-include("parse_eqs.jl")
 
 function __init__()
     @require DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" include("common.jl")
diff --git a/src/common.jl b/src/common.jl
index 3e16d2ecb..053890348 100644
--- a/src/common.jl
+++ b/src/common.jl
@@ -50,20 +50,17 @@ struct TaylorMethodCache{uType, rateType, tTType, uTType} <: OrdinaryDiffEqMutab
     duT::uTType
     uauxT::uTType
     parse_eqs::Ref{Bool}
-    tmpTaylor::uTType
-    arrTaylor::Vector{uTType}
+    rv::TaylorIntegration.RetAlloc
 end
 
-full_cache(c::TaylorMethodCache) = begin
-    tuple(c.u, c.uprev, c.tmp, c.k, c.fsalfirst, c.tT, c.uT, c.duT, c.uauxT, c.parse_eqs,
-        c.tmpTaylor, c.arrTaylor)
-end
+# full_cache(c::TaylorMethodCache) = begin
+#     tuple(c.u, c.uprev, c.tmp, c.k, c.fsalfirst, c.tT, c.uT, c.duT, c.uauxT, c.parse_eqs, c.rv)
+# end
 
 struct TaylorMethodConstantCache{uTType} <: OrdinaryDiffEqConstantCache
     uT::uTType
     parse_eqs::Ref{Bool}
-    tmpTaylor::Vector{uTType}
-    arrTaylor::Vector{Vector{uTType}}
+    rv::TaylorIntegration.RetAlloc{uTType}
 end
 
 function alg_cache(alg::TaylorMethod, u, rate_prototype, uEltypeNoUnits,
@@ -75,8 +72,8 @@ function alg_cache(alg::TaylorMethod, u, rate_prototype, uEltypeNoUnits,
     uT = Taylor1.(u, order)
     duT = zero.(Taylor1.(u, order))
     uauxT = similar(uT)
-    parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(alg.parse_eqs, f, tT, uT, duT, p)
-    TaylorMethodCache(
+    parse_eqs, rv = _determine_parsing!(alg.parse_eqs, f, tT, uT, duT, p)
+    return TaylorMethodCache(
         u,
         uprev,
         similar(u),
@@ -87,8 +84,7 @@ function alg_cache(alg::TaylorMethod, u, rate_prototype, uEltypeNoUnits,
         duT,
         uauxT,
         Ref(parse_eqs),
-        tmpTaylor,
-        arrTaylor
+        rv
         )
 end
 
@@ -103,11 +99,8 @@ function alg_cache(alg::TaylorMethod, u::ArrayPartition, rate_prototype, uEltype
     uT = Taylor1.(u, order)
     duT = zero.(Taylor1.(u, order))
     uauxT = similar(uT)
-    parse_eqs, _, _ = _determine_parsing!(alg.parse_eqs, f, tT, uT, duT, p)
-    ## `tmpT1` must have the same type than `uT` and `arrT1` is a `[tmpT1]`
-    tmpT1 = similar(uT)
-    arrT1 = similar([tmpT1])
-    TaylorMethodCache(
+    parse_eqs, rv = _determine_parsing!(alg.parse_eqs, f, tT, uT, duT, p)
+    return TaylorMethodCache(
         u,
         uprev,
         similar(u),
@@ -118,8 +111,7 @@ function alg_cache(alg::TaylorMethod, u::ArrayPartition, rate_prototype, uEltype
         duT,
         uauxT,
         Ref(parse_eqs),
-        tmpT1,
-        arrT1
+        rv
         )
 end
 
@@ -130,8 +122,8 @@ function alg_cache(alg::TaylorMethod, u, rate_prototype, uEltypeNoUnits,
     tT = Taylor1(typeof(t), order)
     tT[0] = t
     uT = Taylor1(u, order)
-    parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(alg.parse_eqs, f, tT, uT, p)
-    TaylorMethodConstantCache(Taylor1(u, alg.order), Ref(parse_eqs), tmpTaylor, arrTaylor)
+    parse_eqs, rv = _determine_parsing!(alg.parse_eqs, f, tT, uT, p)
+    return TaylorMethodConstantCache(Taylor1(u, alg.order), Ref(parse_eqs), rv)
 end
 
 function initialize!(integrator, c::TaylorMethodConstantCache)
@@ -139,7 +131,7 @@ function initialize!(integrator, c::TaylorMethodConstantCache)
     tT = Taylor1(typeof(t), integrator.alg.order)
     tT[0] = t
     c.uT .= Taylor1(u, tT.order)
-    __jetcoeffs!(Val(c.parse_eqs.x), f, tT, c.uT, p, c.tmpTaylor, c.arrTaylor)
+    __jetcoeffs!(Val(c.parse_eqs.x), f, tT, c.uT, p, c.rv)
     # FSAL stuff
     integrator.kshortsize = 2
     integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
@@ -157,7 +149,7 @@ function perform_step!(integrator,cache::TaylorMethodConstantCache)
     tT[0] = t+dt
     u = evaluate(cache.uT, dt)
     cache.uT[0] = u
-    __jetcoeffs!(Val(cache.parse_eqs.x), f, tT, cache.uT, p, cache.tmpTaylor, cache.arrTaylor)
+    __jetcoeffs!(Val(cache.parse_eqs.x), f, tT, cache.uT, p, cache.rv)
     k = f(u, p, t+dt) # For the interpolation, needs k at the updated point
     integrator.destats.nf += 1
     integrator.fsallast = k
@@ -168,8 +160,8 @@ end
 
 function initialize!(integrator, cache::TaylorMethodCache)
     @unpack u, t, f, p = integrator
-    @unpack k, fsalfirst, tT, uT, duT, uauxT, parse_eqs = cache
-    __jetcoeffs!(Val(parse_eqs.x), f, tT, uT, duT, uauxT, p, cache.tmpTaylor, cache.arrTaylor)
+    @unpack k, fsalfirst, tT, uT, duT, uauxT, parse_eqs, rv = cache
+    __jetcoeffs!(Val(parse_eqs.x), f, tT, uT, duT, uauxT, p, rv)
     # FSAL for interpolation
     integrator.fsalfirst = fsalfirst
     integrator.fsallast = k
@@ -183,14 +175,14 @@ end
 
 function perform_step!(integrator, cache::TaylorMethodCache)
     @unpack t, dt, u, f, p = integrator
-    @unpack k, tT, uT, duT, uauxT, parse_eqs, tmpTaylor, arrTaylor = cache
+    @unpack k, tT, uT, duT, uauxT, parse_eqs, rv = cache
     evaluate!(uT, dt, u)
     tT[0] = t+dt
     for i in eachindex(u)
         @inbounds uT[i][0] = u[i]
         duT[i].coeffs .= zero(duT[i][0])
     end
-    __jetcoeffs!(Val(parse_eqs.x), f, tT, uT, duT, uauxT, p, tmpTaylor, arrTaylor)
+    __jetcoeffs!(Val(parse_eqs.x), f, tT, uT, duT, uauxT, p, rv)
     k = constant_term.(duT) # For the interpolation, needs k at the updated point
     integrator.destats.nf += 1
 end
@@ -206,7 +198,7 @@ function DiffEqBase.solve(
         kwargs...) where
         {uType, tupType, isinplace}
 
-    SciMLBase.unwrapped_f(prob.f)
+    # SciMLBase.unwrapped_f(prob.f)
 
     if verbose
         warned = !isempty(kwargs) && check_keywords(alg, kwargs, warnlist)
@@ -248,13 +240,13 @@ end
 
 # used in continuous callbacks and related methods to update Taylor expansions cache
 function update_jetcoeffs_cache!(u,f,p,cache::TaylorMethodCache)
-    @unpack tT, uT, duT, uauxT, parse_eqs, tmpTaylor, arrTaylor = cache
+    @unpack tT, uT, duT, uauxT, parse_eqs, rv = cache
     @inbounds for i in eachindex(u)
         uT[i][0] = u[i]
         # duT[i].coeffs .= zero(duT[i][0])
         duT[i][0] = zero(uT[i][0])
     end
-    __jetcoeffs!(Val(parse_eqs.x), f, tT, uT, duT, uauxT, p, tmpTaylor, arrTaylor)
+    __jetcoeffs!(Val(parse_eqs.x), f, tT, uT, duT, uauxT, p, rv)
     return nothing
 end
 
@@ -281,22 +273,24 @@ function _ode_addsteps!(k, t, uprev, u, dt, f, p, cache::TaylorMethodCache,
 end
 
 @inline __jetcoeffs!(::Val{false}, f::ODEFunction, t, x::Taylor1{U}, params,
-    tmpTaylor, arrTaylor) where {U} = __jetcoeffs!(Val(false), f.f, t, x, params)
+    rv::RetAlloc{Taylor1{U}}) where {U} = __jetcoeffs!(Val(false), f.f, t, x, params)
 @inline __jetcoeffs!(::Val{true},  f::ODEFunction, t, x::Taylor1{U}, params,
-    tmpTaylor, arrTaylor) where {U} = __jetcoeffs!(Val(true), f.f, t, x, params, tmpTaylor, arrTaylor)
-@inline __jetcoeffs!(::Val{false}, f::ODEFunction, t, x::Array{Taylor1{U},1}, dx, xaux, params,
-    tmpTaylor, arrTaylor) where {U} = __jetcoeffs!(Val(false), f.f, t, x, dx, xaux, params)
-@inline __jetcoeffs!(::Val{true},  f::ODEFunction, t, x::Array{Taylor1{U},1}, dx, xaux, params,
-    tmpTaylor, arrTaylor) where {U} = __jetcoeffs!(Val(true), f.f, t, x, dx, params, tmpTaylor, arrTaylor)
-#
-@inline __jetcoeffs!(::Val{false}, f::DynamicalODEFunction, t, x::Taylor1{U}, params,
-    tmpTaylor, arrTaylor) where {U} = __jetcoeffs!(Val(false), f, t, x, params)
-@inline __jetcoeffs!(::Val{true},  f::DynamicalODEFunction, t, x::Taylor1{U}, params,
-    tmpTaylor, arrTaylor) where {U} = __jetcoeffs!(Val(true), f, t, x, params, tmpTaylor, arrTaylor)
+    rv::RetAlloc{Taylor1{U}}) where {U} = __jetcoeffs!(Val(true), f.f, t, x, params, rv)
+@inline __jetcoeffs!(::Val{false}, f::ODEFunction, t, x::AbstractArray{Taylor1{U},N}, dx, xaux, params,
+    rv::RetAlloc{Taylor1{U}}) where {U,N} = __jetcoeffs!(Val(false), f.f, t, x, dx, xaux, params)
+@inline __jetcoeffs!(::Val{true},  f::ODEFunction, t, x::AbstractArray{Taylor1{U},N}, dx, xaux, params,
+    rv::RetAlloc{Taylor1{U}}) where {U,N} = __jetcoeffs!(Val(true), f.f, t, x, dx, params, rv)
+
+# NOTE: DynamicalODEFunction assumes x is a vector
+# @inline __jetcoeffs!(::Val{false}, f::DynamicalODEFunction, t, x::Taylor1{U}, params,
+#     rv::RetAlloc{Taylor1{U}}) where {U} = __jetcoeffs!(Val(false), f, t, x, params)
+# @inline __jetcoeffs!(::Val{true},  f::DynamicalODEFunction, t, x::Taylor1{U}, params,
+#     rv::RetAlloc{Taylor1{U}}) where {U} = __jetcoeffs!(Val(true), f, t, x, params, rv)
 @inline __jetcoeffs!(::Val{false}, f::DynamicalODEFunction, t, x::ArrayPartition, dx, xaux, params,
-    tmpTaylor, arrTaylor) = jetcoeffs!(f, t, vec(x), vec(dx), xaux, params)
-@inline __jetcoeffs!(::Val{true},  f::DynamicalODEFunction, t, x::ArrayPartition, dx, xaux, params,
-    tmpTaylor, arrTaylor) = __jetcoeffs!(Val(true), f, t, vec(x), vec(dx), params, tmpTaylor, arrTaylor)
+    rv::RetAlloc) = jetcoeffs!(f, t, vec(x), vec(dx), xaux, params)
+# NOTE: `parse_eqs=true` not implemented for `DynamicalODEFunction`
+# @inline __jetcoeffs!(::Val{true},  f::DynamicalODEFunction, t, x::ArrayPartition, dx, xaux, params,
+#     rv::RetAlloc) = __jetcoeffs!(Val(true), f, t, vec(x), vec(dx), params, rv)
 
 _determine_parsing!(parse_eqs::Bool, f::ODEFunction, t, x, params) =
     _determine_parsing!(parse_eqs, f.f, t, x, params)
diff --git a/src/explicitode.jl b/src/explicitode.jl
index 9607a3df0..6c2d64c5f 100644
--- a/src/explicitode.jl
+++ b/src/explicitode.jl
@@ -26,12 +26,12 @@ function jetcoeffs!(eqsdiff::Function, t::Taylor1{T}, x::Taylor1{U}, params) whe
     for ord in 0:order-1
         ordnext = ord+1
 
-        # Set `taux`, `xaux`, auxiliary Taylor1 variables to order `ord`
-        @inbounds taux = Taylor1( t.coeffs[1:ordnext] )
+        # # Set `taux`, `xaux`, auxiliary Taylor1 variables to order `ord`
+        # @inbounds taux = Taylor1( t.coeffs[1:ordnext] )
         @inbounds xaux = Taylor1( x.coeffs[1:ordnext] )
 
         # Equations of motion
-        dx = eqsdiff(xaux, params, taux)
+        dx = eqsdiff(xaux, params, t)
 
         # Recursion relation
         @inbounds x[ordnext] = dx[ord]/ordnext
@@ -67,15 +67,15 @@ function jetcoeffs!(eqsdiff!::Function, t::Taylor1{T},
     for ord in 0:order-1
         ordnext = ord+1
 
-        # Set `taux`, auxiliary Taylor1 variable to order `ord`
-        @inbounds taux = Taylor1( t.coeffs[1:ordnext] )
+        # # Set `taux`, auxiliary Taylor1 variable to order `ord`
+        # @inbounds taux = Taylor1( t.coeffs[1:ordnext] )
         # Set `xaux`, auxiliary vector of Taylor1 to order `ord`
         for j in eachindex(x)
             @inbounds xaux[j] = Taylor1( x[j].coeffs[1:ordnext] )
         end
 
         # Equations of motion
-        eqsdiff!(dx, xaux, params, taux)
+        eqsdiff!(dx, xaux, params, t)
 
         # Recursion relations
         for j in eachindex(x)
@@ -88,29 +88,29 @@ end
 
 # __jetcoeffs
 """
-    __jetcoeffs!(::Val{bool}, f, t, x, params, tmpTaylor, arrTaylor)
-    __jetcoeffs!(::Val{bool}, f, t, x, dx, xaux, params, tmpTaylor, arrTaylor)
+    __jetcoeffs!(::Val{false}, f, t, x, params)
+    __jetcoeffs!(::Val{true}, f, t, x, params, rv)
+    __jetcoeffs!(::Val{false}, f, t, x, dx, xaux, params)
+    __jetcoeffs!(::Val{true}, f, t, x, dx, params, rv)
 
 Chooses a method of [`jetcoeffs!`](@ref) (hard-coded) or the generated by
 [`@taylorize`](@ref)) depending on `Val{bool}` (`bool::Bool`).
 """
 @inline __jetcoeffs!(::Val{false}, f, t, x::Taylor1{U}, params) where {U} =
     jetcoeffs!(f, t, x, params)
-@inline __jetcoeffs!(::Val{true},  f, t, x::Taylor1{U}, params, tmpTaylor, arrTaylor) where {U} =
-    jetcoeffs!(Val(f), t, x, params, tmpTaylor, arrTaylor)
+@inline __jetcoeffs!(::Val{true},  f, t, x::Taylor1{U}, params, rv::RetAlloc{Taylor1{U}}) where {U} =
+    jetcoeffs!(Val(f), t, x, params, rv)
 
 @inline __jetcoeffs!(::Val{false}, f, t, x, dx, xaux, params) =
     jetcoeffs!(f, t, x, dx, xaux, params)
-@inline __jetcoeffs!(::Val{true},  f, t, x, dx, params, tmpTaylor, arrTaylor) =
-    jetcoeffs!(Val(f), t, x, dx, params, tmpTaylor, arrTaylor)
+@inline __jetcoeffs!(::Val{true},  f, t, x, dx, params, rv) =
+    jetcoeffs!(Val(f), t, x, dx, params, rv)
 
 
 # _allocate_jetcoeffs! fallbacks; needed to define the corresponding parsed method
-_allocate_jetcoeffs!(::Taylor1{T}, x::Taylor1{S}, params) where {T,S} =
-    (typeof(x)[], Vector{typeof(x)}[])
-_allocate_jetcoeffs!(::Taylor1{T}, x::AbstractArray{Taylor1{S}, N},
-        ::AbstractArray{Taylor1{S}, N}, params) where {T,S,N} =
-    (eltype(x)[], Vector{eltype(x)}[])
+@inline _allocate_jetcoeffs!(::Taylor1{T}, x::Taylor1{S}, params) where {T,S} = RetAlloc{Taylor1{S}}()
+@inline _allocate_jetcoeffs!(::Taylor1{T}, x::AbstractArray{Taylor1{S}, N},
+        ::AbstractArray{Taylor1{S}, N}, params) where {T,S,N} = RetAlloc{Taylor1{S}}()
 
 
 # stepsize
@@ -238,10 +238,10 @@ function taylorstep!(f, t::Taylor1{T}, x::Taylor1{U}, abstol::T, params) where {
 end
 
 function taylorstep!(f, t::Taylor1{T}, x::Taylor1{U}, abstol::T, params,
-        tmpTaylor::Array{Taylor1{U},1}, arrTaylor) where {T<:Real, U<:Number}
+        rv::RetAlloc{Taylor1{U}}) where {T<:Real, U<:Number}
 
     # Compute the Taylor coefficients
-    __jetcoeffs!(Val(true), f, t, x, params, tmpTaylor, arrTaylor)
+    __jetcoeffs!(Val(true), f, t, x, params, rv)
 
     # Compute the step-size of the integration using `abstol`
     δt = stepsize(x, abstol)
@@ -267,10 +267,10 @@ end
 
 function taylorstep!(f!, t::Taylor1{T}, x::Vector{Taylor1{U}},
         dx::Vector{Taylor1{U}}, abstol::T, params,
-        tmpTaylor::Array{Taylor1{U},1}, arrTaylor) where {T<:Real, U<:Number}
+        rv::RetAlloc{Taylor1{U}}) where {T<:Real, U<:Number}
 
     # Compute the Taylor coefficients
-    __jetcoeffs!(Val(true), f!, t, x, dx, params, tmpTaylor, arrTaylor)
+    __jetcoeffs!(Val(true), f!, t, x, dx, params, rv)
 
     # Compute the step-size of the integration using `abstol`
     δt = stepsize(x, abstol)
@@ -291,12 +291,13 @@ function _determine_parsing!(parse_eqs::Bool, f, t, x, params)
     parse_eqs = parse_eqs && !isempty(methodswith(Val{f}, jetcoeffs!)) &&
                              !isempty(methodswith(Val{f}, _allocate_jetcoeffs!))
 
-    tmpTaylor, arrTaylor = _allocate_jetcoeffs!(t, x, params)
+    # tmpTaylor, arrTaylor = _allocate_jetcoeffs!(t, x, params)
+    rv = _allocate_jetcoeffs!(t, x, params)
 
     if parse_eqs
         try
-            tmpTaylor, arrTaylor = _allocate_jetcoeffs!(Val(f), t, x, params)
-            jetcoeffs!(Val(f), t, x, params, tmpTaylor, arrTaylor)
+            rv = _allocate_jetcoeffs!(Val(f), t, x, params)
+            jetcoeffs!(Val(f), t, x, params, rv)
         catch
             @warn("""Unable to use the parsed method of `jetcoeffs!` for `$f`,
             despite of having `parse_eqs=true`, due to some internal error.
@@ -305,19 +306,19 @@ function _determine_parsing!(parse_eqs::Bool, f, t, x, params)
         end
     end
 
-    return parse_eqs, tmpTaylor, arrTaylor
+    return parse_eqs, rv
 end
 function _determine_parsing!(parse_eqs::Bool, f, t, x, dx, params)
 
     parse_eqs = parse_eqs && !isempty(methodswith(Val{f}, jetcoeffs!)) &&
                              !isempty(methodswith(Val{f}, _allocate_jetcoeffs!))
 
-    tmpTaylor, arrTaylor = _allocate_jetcoeffs!(t, x, dx, params)
+    rv = _allocate_jetcoeffs!(t, x, dx, params)
 
     if parse_eqs
         try
-            tmpTaylor, arrTaylor = _allocate_jetcoeffs!(Val(f), t, x, dx, params)
-            jetcoeffs!(Val(f), t, x, dx, params, tmpTaylor, arrTaylor)
+            rv = _allocate_jetcoeffs!(Val(f), t, x, dx, params)
+            jetcoeffs!(Val(f), t, x, dx, params, rv)
         catch
             @warn("""Unable to use the parsed method of `jetcoeffs!` for `$f`,
             despite of having `parse_eqs=true`, due to some internal error.
@@ -326,7 +327,7 @@ function _determine_parsing!(parse_eqs::Bool, f, t, x, dx, params)
         end
     end
 
-    return parse_eqs, tmpTaylor, arrTaylor
+    return parse_eqs, rv
 end
 
 
@@ -341,13 +342,13 @@ for V in (:(Val{true}), :(Val{false}))
             x = Taylor1( x0, order )
 
             # Determine if specialized jetcoeffs! method exists
-            parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(parse_eqs, f, t, x, params)
+            parse_eqs, rv = _determine_parsing!(parse_eqs, f, t, x, params)
 
             if parse_eqs
                 # Re-initialize the Taylor1 expansions
                 t = t0 + Taylor1( T, order )
                 x = Taylor1( x0, order )
-                return _taylorinteg!(f, t, x, x0, t0, tmax, abstol, tmpTaylor, arrTaylor,
+                return _taylorinteg!(f, t, x, x0, t0, tmax, abstol, rv,
                     $V(), params, maxsteps=maxsteps)
             else
                 return _taylorinteg!(f, t, x, x0, t0, tmax, abstol,
@@ -406,7 +407,7 @@ for V in (:(Val{true}), :(Val{false}))
             end
         end
         function _taylorinteg!(f, t::Taylor1{T}, x::Taylor1{U},
-                x0::U, t0::T, tmax::T, abstol::T, tmpTaylor, arrTaylor, ::$V, params;
+                x0::U, t0::T, tmax::T, abstol::T, rv::RetAlloc{Taylor1{U}}, ::$V, params;
                 maxsteps::Int=500) where {T<:Real, U<:Number}
 
             # Allocation
@@ -428,7 +429,7 @@ for V in (:(Val{true}), :(Val{false}))
 
             # Integration
             while sign_tstep*t0 < sign_tstep*tmax
-                δt = taylorstep!(f, t, x, abstol, params, tmpTaylor, arrTaylor) # δt is positive!
+                δt = taylorstep!(f, t, x, abstol, params, rv) # δt is positive!
                 # Below, δt has the proper sign according to the direction of the integration
                 δt = sign_tstep * min(δt, sign_tstep*(tmax-t0))
                 x0 = evaluate(x, δt) # new initial condition
@@ -471,13 +472,13 @@ for V in (:(Val{true}), :(Val{false}))
             end
 
             # Determine if specialized jetcoeffs! method exists
-            parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(parse_eqs, f!, t, x, dx, params)
+            parse_eqs, rv = _determine_parsing!(parse_eqs, f!, t, x, dx, params)
 
             if parse_eqs
                 # Re-initialize the Taylor1 expansions
                 t = t0 + Taylor1( T, order )
                 x .= Taylor1.( q0, order )
-                return _taylorinteg!(f!, t, x, dx, q0, t0, tmax, abstol, tmpTaylor, arrTaylor,
+                return _taylorinteg!(f!, t, x, dx, q0, t0, tmax, abstol, rv,
                     $V(), params, maxsteps=maxsteps)
             else
                 return _taylorinteg!(f!, t, x, dx, q0, t0, tmax, abstol,
@@ -548,7 +549,7 @@ for V in (:(Val{true}), :(Val{false}))
         end
 
         function _taylorinteg!(f!, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{Taylor1{U},1},
-                q0::Array{U,1}, t0::T, tmax::T, abstol::T, tmpTaylor, arrTaylor, ::$V, params;
+                q0::Array{U,1}, t0::T, tmax::T, abstol::T, rv::RetAlloc{Taylor1{U}}, ::$V, params;
                 maxsteps::Int=500) where {T<:Real, U<:Number}
 
             # Initialize the vector of Taylor1 expansions
@@ -574,7 +575,7 @@ for V in (:(Val{true}), :(Val{false}))
             # Integration
             nsteps = 1
             while sign_tstep*t0 < sign_tstep*tmax
-                δt = taylorstep!(f!, t, x, dx, abstol, params, tmpTaylor, arrTaylor) # δt is positive!
+                δt = taylorstep!(f!, t, x, dx, abstol, params, rv) # δt is positive!
                 # Below, δt has the proper sign according to the direction of the integration
                 δt = sign_tstep * min(δt, sign_tstep*(tmax-t0))
                 evaluate!(x, δt, x0) # new initial condition
@@ -713,10 +714,10 @@ function taylorinteg(f, x0::U, trange::AbstractVector{T},
     x = Taylor1( x0, order )
 
     # Determine if specialized jetcoeffs! method exists
-    parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(parse_eqs, f, t, x, params)
+    parse_eqs, rv = _determine_parsing!(parse_eqs, f, t, x, params)
 
     if parse_eqs
-        return _taylorinteg!(f, t, x, x0, trange, abstol, tmpTaylor, arrTaylor,
+        return _taylorinteg!(f, t, x, x0, trange, abstol, rv,
             params, maxsteps=maxsteps)
     else
         return _taylorinteg!(f, t, x, x0, trange, abstol, params, maxsteps=maxsteps)
@@ -771,7 +772,7 @@ function _taylorinteg!(f, t::Taylor1{T}, x::Taylor1{U}, x0::U, trange::AbstractV
     return xv
 end
 function _taylorinteg!(f, t::Taylor1{T}, x::Taylor1{U}, x0::U, trange::AbstractVector{T},
-        abstol::T, tmpTaylor, arrTaylor, params; maxsteps::Int=500) where {T<:Real, U<:Number}
+        abstol::T, rv::RetAlloc{Taylor1{U}}, params; maxsteps::Int=500) where {T<:Real, U<:Number}
 
     # Allocation
     nn = length(trange)
@@ -788,7 +789,7 @@ function _taylorinteg!(f, t::Taylor1{T}, x::Taylor1{U}, x0::U, trange::AbstractV
     iter = 2
     nsteps = 1
     while sign_tstep*t0 < sign_tstep*tmax
-        δt = taylorstep!(f, t, x, abstol, params, tmpTaylor, arrTaylor)# δt is positive!
+        δt = taylorstep!(f, t, x, abstol, params, rv)# δt is positive!
         # Below, δt has the proper sign according to the direction of the integration
         δt = sign_tstep * min(δt, sign_tstep*(tmax-t0))
         x0 = evaluate(x, δt) # new initial condition
@@ -837,10 +838,10 @@ function taylorinteg(f!, q0::Array{U,1}, trange::AbstractVector{T},
     end
 
     # Determine if specialized jetcoeffs! method exists
-    parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(parse_eqs, f!, t, x, dx, params)
+    parse_eqs, rv = _determine_parsing!(parse_eqs, f!, t, x, dx, params)
 
     if parse_eqs
-        return _taylorinteg!(f!, t, x, dx, q0, trange, abstol, tmpTaylor, arrTaylor,
+        return _taylorinteg!(f!, t, x, dx, q0, trange, abstol, rv,
             params, maxsteps=maxsteps)
     else
         return _taylorinteg!(f!, t, x, dx, q0, trange, abstol, params, maxsteps=maxsteps)
@@ -909,7 +910,7 @@ function _taylorinteg!(f!, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{Tayl
     return transpose(xv)
 end
 function _taylorinteg!(f!, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{Taylor1{U},1},
-        q0::Array{U,1}, trange::AbstractVector{T}, abstol::T, tmpTaylor, arrTaylor, params;
+        q0::Array{U,1}, trange::AbstractVector{T}, abstol::T, rv::RetAlloc{Taylor1{U}}, params;
         maxsteps::Int=500) where {T<:Real, U<:Number}
 
     # Allocation
@@ -935,7 +936,7 @@ function _taylorinteg!(f!, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{Tayl
     iter = 2
     nsteps = 1
     while sign_tstep*t0 < sign_tstep*tmax
-        δt = taylorstep!(f!, t, x, dx, abstol, params, tmpTaylor, arrTaylor) # δt is positive!
+        δt = taylorstep!(f!, t, x, dx, abstol, params, rv) # δt is positive!
         # Below, δt has the proper sign according to the direction of the integration
         δt = sign_tstep * min(δt, sign_tstep*(tmax-t0))
         evaluate!(x, δt, x0) # new initial condition
diff --git a/src/lyapunovspectrum.jl b/src/lyapunovspectrum.jl
index b9879bff4..bd50f4307 100644
--- a/src/lyapunovspectrum.jl
+++ b/src/lyapunovspectrum.jl
@@ -190,11 +190,10 @@ function lyap_taylorstep!(f!, t::Taylor1{T}, x::Vector{Taylor1{U}},
 
     return δt
 end
-function lyap_taylorstep!(f!, t::Taylor1{T}, x::Vector{Taylor1{U}},
-    dx::Vector{Taylor1{U}}, #xaux::Vector{Taylor1{U}},
+function lyap_taylorstep!(f!, t::Taylor1{T}, x::Vector{Taylor1{U}}, dx::Vector{Taylor1{U}},
     δx::Array{TaylorN{Taylor1{U}},1}, dδx::Array{TaylorN{Taylor1{U}},1},
     jac::Array{Taylor1{U},2}, abstol::T, _δv::Vector{TaylorN{Taylor1{U}}},
-    varsaux::Array{Taylor1{U},3}, params, tmpTaylor, arrTaylor,
+    varsaux::Array{Taylor1{U},3}, params, rv::RetAlloc{Taylor1{U}},
     jacobianfunc! =nothing) where {T<:Real, U<:Number}
 
     # Dimensions of phase-space: dof
@@ -202,7 +201,7 @@ function lyap_taylorstep!(f!, t::Taylor1{T}, x::Vector{Taylor1{U}},
     dof = length(δx)
 
     # Compute the Taylor coefficients associated to trajectory
-    __jetcoeffs!(Val(true), f!, t, view(x, 1:dof), view(dx, 1:dof), params, tmpTaylor, arrTaylor)
+    __jetcoeffs!(Val(true), f!, t, view(x, 1:dof), view(dx, 1:dof), params, rv)
 
     # Compute stability matrix
     stabilitymatrix!(f!, t, x, δx, dδx, jac, _δv, params, jacobianfunc!)
@@ -258,14 +257,14 @@ function lyap_taylorinteg(f!, q0::Array{U,1}, t0::T, tmax::T,
     end
 
     # Determine if specialized jetcoeffs! method exists
-    parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(parse_eqs, f!, t,
-                                        view(x, 1:dof), view(dx, 1:dof), params)
+    parse_eqs, rv = _determine_parsing!(parse_eqs, f!, t,
+                    view(x, 1:dof), view(dx, 1:dof), params)
 
     if parse_eqs
         # Re-initialize the Taylor1 expansions
         t = t0 + Taylor1( T, order )
         x .= Taylor1.( x0, order )
-        return _lyap_taylorinteg!(f!, t, x, dx, q0, t0, tmax, abstol, jt, _dv, tmpTaylor, arrTaylor,
+        return _lyap_taylorinteg!(f!, t, x, dx, q0, t0, tmax, abstol, jt, _dv, rv,
             params, jacobianfunc!, maxsteps=maxsteps)
     else
         return _lyap_taylorinteg!(f!, t, x, dx, q0, t0, tmax, abstol, jt, _dv,
@@ -352,8 +351,7 @@ end
 
 function _lyap_taylorinteg!(f!, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{Taylor1{U},1},
         q0::Array{U,1}, t0::T, tmax::T, abstol::T, jt::Matrix{U}, _δv::Array{TaylorN{Taylor1{U}},1},
-        tmpTaylor::Array{Taylor1{U},1}, arrTaylor,
-        params, jacobianfunc!; maxsteps::Int=500) where {T<:Real, U<:Number}
+        rv::RetAlloc{Taylor1{U}}, params, jacobianfunc!; maxsteps::Int=500) where {T<:Real, U<:Number}
 
     # Allocation
     order = get_order(t)
@@ -392,7 +390,7 @@ function _lyap_taylorinteg!(f!, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array
     nsteps = 1
     while sign_tstep*t0 < sign_tstep*tmax
         δt = lyap_taylorstep!(f!, t, x, dx, δx, dδx, jac, abstol, _δv,
-            varsaux, params, tmpTaylor, arrTaylor, jacobianfunc!) # δt is positive!
+            varsaux, params, rv, jacobianfunc!) # δt is positive!
         # Below, δt has the proper sign according to the direction of the integration
         δt = sign_tstep * min(δt, sign_tstep*(tmax-t0))
         evaluate!(x, δt, x0) # Update x0
@@ -455,14 +453,15 @@ function lyap_taylorinteg(f!, q0::Array{U,1}, trange::AbstractVector{T},
     end
 
     # Determine if specialized jetcoeffs! method exists
-    parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(parse_eqs, f!, t,
-                                        view(x, 1:dof), view(dx, 1:dof), params)
+    parse_eqs, rv = _determine_parsing!(parse_eqs, f!, t,
+                    view(x, 1:dof), view(dx, 1:dof), params)
 
     if parse_eqs
         # Re-initialize the Taylor1 expansions
         t = trange[1] + Taylor1( T, order )
         x .= Taylor1.( x0, order )
-        return _lyap_taylorinteg!(f!, t, x, dx, q0, trange, abstol, jt, _dv, tmpTaylor, arrTaylor,
+        tmpTaylor, arrTaylor = rv.v0, rv.v1
+        return _lyap_taylorinteg!(f!, t, x, dx, q0, trange, abstol, jt, _dv, rv,
             params, jacobianfunc!, maxsteps=maxsteps)
     else
         return _lyap_taylorinteg!(f!, t, x, dx, q0, trange, abstol, jt, _dv,
@@ -579,8 +578,7 @@ end
 
 function _lyap_taylorinteg!(f!, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{Taylor1{U},1},
         q0::Array{U,1}, trange::AbstractVector{T}, abstol::T, jt::Matrix{U}, _δv::Array{TaylorN{Taylor1{U}},1},
-        tmpTaylor::Vector{Taylor1{U}}, arrTaylor,
-        params, jacobianfunc!; maxsteps::Int=500) where {T<:Real, U<:Number}
+        rv::RetAlloc{Taylor1{U}}, params, jacobianfunc!; maxsteps::Int=500) where {T<:Real, U<:Number}
 
     # Allocation
     order = get_order(t)
@@ -623,7 +621,7 @@ function _lyap_taylorinteg!(f!, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array
     nsteps = 1
     while sign_tstep*t0 < sign_tstep*tmax
         δt = lyap_taylorstep!(f!, t, x, dx, δx, dδx, jac, abstol, _δv,
-            varsaux, params, tmpTaylor, arrTaylor, jacobianfunc!) # δt is positive!
+            varsaux, params, rv, jacobianfunc!) # δt is positive!
         # Below, δt has the proper sign according to the direction of the integration
         δt = sign_tstep * min(δt, sign_tstep*(tmax-t0))
         evaluate!(x, δt, x0) # Update x0
diff --git a/src/parse_eqs.jl b/src/parse_eqs.jl
index a5a6b1d7a..eb319314c 100644
--- a/src/parse_eqs.jl
+++ b/src/parse_eqs.jl
@@ -15,6 +15,10 @@ Mutable struct that contains all the bookkeeping vectors/dictionaries used withi
     - `d_decl`     : Dictionary declared arrays
     - `v_newvars`  : Symbols of auxiliary indexed vars
     - `v_arraydecl`: Symbols which are explicitly declared as Array or Vector
+    - `v_array1`   : Symbols which are explicitly declared as Array{Taylor1{T},1}
+    - `v_array2`   : Symbols which are explicitly declared as Array{Taylor1{T},2}
+    - `v_array3`   : Symbols which are explicitly declared as Array{Taylor1{T},3}
+    - `v_array4`   : Symbols which are explicitly declared as Array{Taylor1{T},4}
     - `v_preamb`   : Symbols or Expr used in the preamble (declarations, etc)
     - `retvar`     : *Guessed* returned variable, which defines the LHS of the ODEs
 
@@ -25,15 +29,50 @@ mutable struct BookKeeping
     d_decl   :: Dict{Symbol, Expr}
     v_newvars   :: Vector{Symbol}
     v_arraydecl :: Vector{Symbol}
+    v_array1 :: Vector{Symbol}
+    v_array2 :: Vector{Symbol}
+    v_array3 :: Vector{Symbol}
+    v_array4 :: Vector{Symbol}
     v_preamb :: Vector{Union{Symbol,Expr}}
     retvar   :: Symbol
 
     function BookKeeping()
         return new(Dict{Symbol, Expr}(), Dict{Union{Symbol,Expr}, Number}(),
-            Dict{Symbol, Expr}(), Symbol[], Symbol[], Union{Symbol,Expr}[], :nothing)
+            Dict{Symbol, Expr}(), Symbol[], Symbol[], Symbol[], Symbol[], Symbol[], Symbol[],
+            Union{Symbol,Expr}[], :nothing)
     end
 end
 
+
+"""
+    RetAlloc{Taylor1{T}}
+
+Struct related to the returned variables that are pre-allocated when
+`@taylorize` is used.
+    - `v0`   : Vector{Taylor1{T}}
+    - `v1`   : Vector{Vector{Taylor1{T}}}
+
+"""
+struct RetAlloc{T <: Number}
+    v0 :: Array{T,1}
+    v1 :: Vector{Array{T,1}}
+    v2 :: Vector{Array{T,2}}
+    v3 :: Vector{Array{T,3}}
+    v4 :: Vector{Array{T,4}}
+
+    function RetAlloc{T}() where {T}
+        v1 = Array{T,1}(undef, 0)
+        return new(v1, [v1], [Array{T,2}(undef, 0, 0)], [Array{T,3}(undef, 0, 0, 0)],
+            [Array{T,4}(undef, 0, 0, 0, 0)])
+    end
+
+    function RetAlloc{T}(v0::Array{T,1}, v1::Vector{Array{T,1}},
+            v2::Vector{Array{T,2}}, v3::Vector{Array{T,3}}, v4::Vector{Array{T,4}}) where {T}
+        return new(v0, v1, v2, v3, v4)
+    end
+end
+
+
 """
 `inbookkeeping(v, bkkeep::BookKeeping)`
 
@@ -49,7 +88,7 @@ Checks if `v` is declared in `bkkeep`, considering the `d_indx`, `v_newvars` and
 # The (irrelevant) `nothing` below is there to have a `:block` Expr; it is deleted later
 const _HEAD_PARSEDFN_SCALAR = sanitize(:(
     function TaylorIntegration.jetcoeffs!(::Val{__fn}, __tT::Taylor1{_T}, __x::Taylor1{_S},
-            __params, __tmpTaylor, __arrTaylor) where {_T<:Real, _S<:Number}
+            __params, __ralloc::TaylorIntegration.RetAlloc{Taylor1{_S}}) where {_T<:Real, _S<:Number}
         order = __tT.order
         nothing
     end)
@@ -58,7 +97,7 @@ const _HEAD_PARSEDFN_SCALAR = sanitize(:(
 const _HEAD_PARSEDFN_VECTOR = sanitize(:(
     function TaylorIntegration.jetcoeffs!(::Val{__fn}, __tT::Taylor1{_T},
             __x::AbstractArray{Taylor1{_S}, _N}, __dx::AbstractArray{Taylor1{_S}, _N},
-            __params, __tmpTaylor, __arrTaylor) where {_T<:Real, _S<:Number, _N}
+            __params, __ralloc::TaylorIntegration.RetAlloc{Taylor1{_S}}) where {_T<:Real, _S<:Number, _N}
         order = __tT.order
         nothing
     end)
@@ -89,7 +128,7 @@ const _DECL_ARRAY = sanitize( Expr(:block,
 
 
 """
-`_make_parsed_jetcoeffs( ex, debug=false )`
+`_make_parsed_jetcoeffs( ex )`
 
 This function constructs the expressions of two new methods, the first equivalent to the
 differential equations (jetcoeffs!), which exploits the mutating functions of TaylorSeries.jl,
@@ -97,61 +136,51 @@ and the second one (_allocate_jetcoeffs) preallocates any auxiliary `Taylor1` or
 `Vector{Taylor1{T}}` needed.
 
 """
-function _make_parsed_jetcoeffs(ex::Expr, debug=false)
+function _make_parsed_jetcoeffs(ex::Expr)
     # Extract the name, args and body of the function
     # `fn` name of the function having the ODEs
     # `fnargs` arguments of the function
     # `fnbody` future transformed function body
     fn, fnargs, fnbody = _extract_parts(ex)
-    debug && (println("****** _extract_parts ******");
-        @show(fn, fnargs, fnbody); println())
 
     # Set up the Expr for the new functions
     new_jetcoeffs, new_allocjetcoeffs = _newhead(fn, fnargs)
 
-    # Expr for the for-loop block for the recursion (of the `x` variable)
-    forloopblock = Expr(:for, :(ord = 1:order-1), Expr(:block, :(ordnext = ord + 1)) )
-
     # Transform the graph representation of the body of the functions:
     # defspreamble: inicializations used for the zeroth order (preamble)
     # defsprealloc: definitions (declarations) of auxiliary Taylor1's
     # fnbody: transformed function body, using mutating functions from TaylorSeries;
     #         used later within the recursion loop
     # bkkeep: book-keeping structure having info of the variables
-    defspreamble, defsprealloc, fnbody, bkkeep = _preamble_body(fnbody, fnargs, debug)
-    debug && (println("****** _preamble_body ******");
-        @show(defspreamble); println(); @show(defsprealloc); println();
-        @show(fnbody); println(); @show(bkkeep); println())
+    defspreamble, defsprealloc, fnbody, bkkeep = _preamble_body(fnbody, fnargs)
 
     # Create body of recursion loop; temporary assignements may be needed.
     # rec_preamb: recursion loop for the preamble (first order correction)
     # rec_fnbody: recursion loop for the body-function (recursion loop for higher orders)
     rec_preamb, rec_fnbody = _recursionloop(fnargs, bkkeep)
-    debug && (println("****** _recursionloop ******");
-        @show(rec_preamb); println(); @show(rec_fnbody); println())
 
+    # Expr for the for-loop block for the recursion (of the `x` variable)
+    forloopblock = Expr(:for, :(ord = 1:order-1), Expr(:block, :(ordnext = ord + 1)) )
     # Add rec_fnbody to forloopblock
     push!(forloopblock.args[2].args, fnbody.args[1].args..., rec_fnbody)
 
-    # Add preamble and recursion body to `new_allocjetcoeffs`
+    # Add preamble and recursion body to `new_jetcoeffs`
     push!(new_jetcoeffs.args[2].args, defspreamble..., rec_preamb)
 
-    # Push preamble and forloopblock to `new_allocjetcoeffs` and return line
+    # Push preamble and forloopblock to `new_jetcoeffs` and return line
     push!(new_jetcoeffs.args[2].args, forloopblock, Meta.parse("return nothing"))
 
+    # Split v_arraydecl according to the number of indices
+    _split_arraydecl!(bkkeep)
+
+    # Add allocated variable definitions to `new_jetcoeffs`, to make it more human readable
+    _allocated_defs!(new_jetcoeffs, bkkeep)
+
     # Define the expressions of the returned vectors in `new_allocjetcoeffs`
     push!(new_allocjetcoeffs.args[2].args, defsprealloc...)
-    if isempty(bkkeep.v_newvars)
-        ret_defsprealloc = :(Taylor1{_S}[])
-    else
-        ret_defsprealloc = :([$(bkkeep.v_newvars...),])
-    end
-    if isempty(bkkeep.v_arraydecl)
-        ret_tmparrays = :([Vector{Taylor1{_S}}(undef, 0),])
-    else
-        ret_tmparrays = :([$(bkkeep.v_arraydecl...),])
-    end
-    ret_ret = :(return $(ret_defsprealloc), $(ret_tmparrays))
+
+    # Define returned expression for `new_allocjetcoeffs`
+    ret_ret = _returned_expr(bkkeep)
 
     # Add return line to `new_allocjetcoeffs`
     push!(new_allocjetcoeffs.args[2].args, ret_ret)
@@ -172,21 +201,10 @@ function _make_parsed_jetcoeffs(ex::Expr, debug=false)
             Dict(:__tT => fnargs[4], :__params => fnargs[3],
                 :(__x) => fnargs[2], :(__dx) => fnargs[1], :(__fn) => fn ))
 
-    else
-        # A priori this is not needed
-        throw(ArgumentError("Wrong number of arguments in `fnargs`"))
-    end
-
-    # Bring back old variable-names to `new_jecoeffs`; includes proper renaming
-    # of the temporal arrays
-    dd = Dict{Symbol, Expr}()
-    @inbounds for (ind, vnew) in enumerate(bkkeep.v_newvars)
-        merge!(dd, Dict(vnew => Expr(:ref, :__tmpTaylor, ind)))
+    # else
+    #     # A priori this is not needed
+    #     throw(ArgumentError("Wrong number of arguments in `fnargs`"))
     end
-    @inbounds for (ind, vnew) in enumerate(bkkeep.v_arraydecl)
-        merge!(dd, Dict(vnew => Expr(:ref, :__arrTaylor, ind)))
-    end
-    new_jetcoeffs = subs(new_jetcoeffs, dd)
 
     return new_jetcoeffs, new_allocjetcoeffs
 end
@@ -309,7 +327,7 @@ end
 
 
 """
-`_preamble_body(fnbody, fnargs, debug=false)`
+`_preamble_body(fnbody, fnargs)`
 
 Returns expressions for the preamble, the declaration of
 arrays, the body and the bookkeeping struct, which will be used to build
@@ -318,21 +336,17 @@ the original function (already adapted), `fnargs` is a vector of symbols
 of the original diferential equations function.
 
 """
-function _preamble_body(fnbody, fnargs, debug=false)
+function _preamble_body(fnbody, fnargs)
     # Inicialize BookKeeping struct
     bkkeep = BookKeeping()
 
     # Rename vars to have the body in non-indexed form; bkkeep has different entries
     # for bookkeeping variables/symbolds, including indexed ones
     fnbody, bkkeep.d_indx = _rename_indexedvars(fnbody)
-    debug && (println("------ _rename_indexedvars ------");
-        @show(fnbody); println(); @show(bkkeep.d_indx); println())
 
     # Create `newfnbody` which corresponds to `fnbody`, cleaned (without irrelevant comments)
     # and with all new variables in place; bkkeep.d_indx is updated
     newfnbody = _newfnbody!(fnbody, fnargs, bkkeep)
-    debug && (println("------ _newfnbody ------");
-        @show(newfnbody); println(); @show(bkkeep); println())
 
     # Parse `newfnbody` and create `prepreamble` and `prealloc`, updating `bkkeep`.
     # These objects use the mutating functions from TaylorSeries.
@@ -347,9 +361,6 @@ function _preamble_body(fnbody, fnargs, debug=false)
     newfnbody = subs(newfnbody, bkkeep.d_assign)
     preamble = subs(preamble, bkkeep.d_assign)
     prealloc = subs(prealloc, bkkeep.d_assign)
-    debug && (println("------ _parse_newfnbody! ------");
-        @show(newfnbody); println(); @show(preamble); println();
-        @show(prealloc); println(); @show(bkkeep); println())
 
     # Include the assignement of indexed auxiliary variables
     defsprealloc = _defs_allocs!(prealloc, fnargs, bkkeep, false)
@@ -361,10 +372,6 @@ function _preamble_body(fnbody, fnargs, debug=false)
     # Define retvar; for scalar eqs is the last entry included in v_newvars
     bkkeep.retvar = length(fnargs) == 3 ? subs(bkkeep.v_newvars[end], bkkeep.d_indx) : fnargs[1]
 
-    debug && (println("------ _defs_allocs! ------");
-        @show(defsprealloc); println(); @show(preamble); println();
-        @show(newfnbody); println(); @show(bkkeep); println())
-
     return defspreamble, defsprealloc, newfnbody, bkkeep
 end
 
@@ -1053,6 +1060,103 @@ function _recursionloop(fnargs, bkkeep::BookKeeping)
     return rec_preamb, rec_fnbody
 end
 
+
+"""
+`_split_arraydecl!(bkkeep)`
+
+Split bkkeep.v_arraydecl in the vector (bkkeep.v_array1), matrix (bkkeep.v_array2), etc,
+to properly construct the `RetAlloc` variable.
+"""
+function _split_arraydecl!(bkkeep::BookKeeping)
+    for s in bkkeep.v_arraydecl
+        for v in values(bkkeep.d_indx)
+            if v.head == :ref && v.args[1] == s
+                if length(v.args) == 2
+                    push!(bkkeep.v_array1, s)
+                    break
+                elseif length(v.args) == 3
+                    push!(bkkeep.v_array2, s)
+                    break
+                elseif length(v.args) == 4
+                    push!(bkkeep.v_array3, s)
+                    break
+                elseif length(v.args) == 5
+                    push!(bkkeep.v_array4, s)
+                    break
+                else
+                    error("Error: `@taylorize` allows only to parse up tp 5-index arrays")
+                end
+            end
+        end
+    end
+    return nothing
+end
+
+
+"""
+`_allocated_defs!(new_jetcoeffs, bkkeep)`
+
+Add allocated variable definitions to `new_jetcoeffs`, to make it more human readable.
+"""
+function _allocated_defs!(new_jetcoeffs::Expr, bkkeep::BookKeeping)
+    tmp_defs = [popfirst!(new_jetcoeffs.args[2].args)]
+    @inbounds for (ind, vnew) in enumerate(bkkeep.v_newvars)
+        push!(tmp_defs, :($(vnew) = __ralloc.v0[$(ind)]))
+    end
+    @inbounds for (ind, vnew) in enumerate(bkkeep.v_array1)
+        push!(tmp_defs, :($(vnew) = __ralloc.v1[$(ind)]))
+    end
+    @inbounds for (ind, vnew) in enumerate(bkkeep.v_array2)
+        push!(tmp_defs, :($(vnew) = __ralloc.v2[$(ind)]))
+    end
+    @inbounds for (ind, vnew) in enumerate(bkkeep.v_array3)
+        push!(tmp_defs, :($(vnew) = __ralloc.v3[$(ind)]))
+    end
+    @inbounds for (ind, vnew) in enumerate(bkkeep.v_array4)
+        push!(tmp_defs, :($(vnew) = __ralloc.v4[$(ind)]))
+    end
+    prepend!(new_jetcoeffs.args[2].args, tmp_defs)
+    return nothing
+end
+
+
+"""
+`_returned_expr(bkkeep)`
+
+Constructs the expression to be returned by `TaylorIntegration._allocate_jetcoeffs!`
+"""
+function _returned_expr(bkkeep::BookKeeping)
+    if isempty(bkkeep.v_newvars)
+        retv0 = :(Taylor1{_S}[])
+    else
+        retv0 = :([$(bkkeep.v_newvars...),])
+    end
+    if isempty(bkkeep.v_array1)
+        retv1 = :([Array{Taylor1{_S},1}(undef, 0),])
+    else
+        retv1 = :([$(bkkeep.v_array1...),])
+    end
+    if isempty(bkkeep.v_array2)
+        retv2 = :([Array{Taylor1{_S},2}(undef, 0, 0),])
+    else
+        retv2 = :([$(bkkeep.v_array2...),])
+    end
+    if isempty(bkkeep.v_array3)
+        retv3 = :([Array{Taylor1{_S},3}(undef, 0, 0, 0),])
+    else
+        retv3 = :([$(bkkeep.v_array3...),])
+    end
+    if isempty(bkkeep.v_array4)
+        retv4 = :([Array{Taylor1{_S},4}(undef, 0, 0, 0, 0),])
+    else
+        retv4 = :([$(bkkeep.v_array4...),])
+    end
+
+    return :(return TaylorIntegration.RetAlloc{Taylor1{_S}}(
+        $(retv0), $(retv1), $(retv2), $(retv3), $(retv4)))
+end
+
+
 """
 `@taylorize expr`
 
diff --git a/src/rootfinding.jl b/src/rootfinding.jl
index 929535af6..592f36769 100644
--- a/src/rootfinding.jl
+++ b/src/rootfinding.jl
@@ -174,13 +174,13 @@ function taylorinteg(f!, g, q0::Array{U,1}, t0::T, tmax::T,
     end
 
     # Determine if specialized jetcoeffs! method exists
-    parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(parse_eqs, f!, t, x, dx, params)
+    parse_eqs, rv = _determine_parsing!(parse_eqs, f!, t, x, dx, params)
 
     if parse_eqs
         # Re-initialize the Taylor1 expansions
         t = t0 + Taylor1( T, order )
         x .= Taylor1.( q0, order )
-        return _taylorinteg!(f!, g, t, x, dx, q0, t0, tmax, abstol, tmpTaylor, arrTaylor, params,
+        return _taylorinteg!(f!, g, t, x, dx, q0, t0, tmax, abstol, rv, params,
             maxsteps=maxsteps, eventorder=eventorder, newtoniter=newtoniter, nrabstol=nrabstol)
     else
         return _taylorinteg!(f!, g, t, x, dx, q0, t0, tmax, abstol,params,
@@ -257,7 +257,7 @@ function _taylorinteg!(f!, g, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{T
     return view(tv,1:nsteps), view(transpose(view(xv,:,1:nsteps)),1:nsteps,:), view(tvS,1:nevents-1), view(transpose(view(xvS,:,1:nevents-1)),1:nevents-1,:), view(gvS,1:nevents-1)
 end
 function _taylorinteg!(f!, g, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{Taylor1{U},1},
-        q0::Array{U,1}, t0::T, tmax::T, abstol::T, tmpTaylor, arrTaylor, params;
+        q0::Array{U,1}, t0::T, tmax::T, abstol::T, rv::RetAlloc{Taylor1{U}}, params;
         maxsteps::Int=500, eventorder::Int=0, newtoniter::Int=10, nrabstol::T=eps(T)) where {T <: Real,U <: Number}
 
     # Allocation
@@ -298,7 +298,7 @@ function _taylorinteg!(f!, g, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{T
     while sign_tstep*t0 < sign_tstep*tmax
         δt_old = δt
         # δt = taylorstep!(f!, t, x, dx, xaux, abstol, params) # δt is positive!
-        δt = taylorstep!(f!, t, x, dx, abstol, params, tmpTaylor, arrTaylor) # δt is positive!
+        δt = taylorstep!(f!, t, x, dx, abstol, params, rv) # δt is positive!
         # Below, δt has the proper sign according to the direction of the integration
         δt = sign_tstep * min(δt, sign_tstep*(tmax-t0))
         evaluate!(x, δt, x0) # new initial condition
@@ -349,13 +349,13 @@ function taylorinteg(f!, g, q0::Array{U,1}, trange::AbstractVector{T},
     end
 
     # Determine if specialized jetcoeffs! method exists
-    parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(parse_eqs, f!, t, x, dx, params)
+    parse_eqs, rv = _determine_parsing!(parse_eqs, f!, t, x, dx, params)
 
     if parse_eqs
         # Re-initialize the Taylor1 expansions
         t = t0 + Taylor1( T, order )
         x .= Taylor1.(q0, order)
-        return _taylorinteg!(f!, g, t, x, dx, q0, trange, abstol, tmpTaylor, arrTaylor, params,
+        return _taylorinteg!(f!, g, t, x, dx, q0, trange, abstol, rv, params,
             maxsteps=maxsteps, eventorder=eventorder, newtoniter=newtoniter, nrabstol=nrabstol)
     else
         return _taylorinteg!(f!, g, t, x, dx, q0, trange, abstol,params,
@@ -446,7 +446,7 @@ function _taylorinteg!(f!, g, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{T
 end
 
 function _taylorinteg!(f!, g, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{Taylor1{U},1},
-        q0::Array{U,1}, trange::AbstractVector{T}, abstol::T, tmpTaylor, arrTaylor, params;
+        q0::Array{U,1}, trange::AbstractVector{T}, abstol::T, rv::RetAlloc{Taylor1{U}}, params;
         maxsteps::Int=500, eventorder::Int=0, newtoniter::Int=10, nrabstol::T=eps(T)) where {T <: Real,U <: Number}
 
     # Allocation
@@ -487,7 +487,7 @@ function _taylorinteg!(f!, g, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{T
     nevents = 1 #number of detected events
     while sign_tstep*t0 < sign_tstep*tmax
         δt_old = δt
-        δt = taylorstep!(f!, t, x, dx, abstol, params, tmpTaylor, arrTaylor) # δt is positive!
+        δt = taylorstep!(f!, t, x, dx, abstol, params, rv) # δt is positive!
         # Below, δt has the proper sign according to the direction of the integration
         δt = sign_tstep * min(δt, sign_tstep*(tmax-t0))
         evaluate!(x, δt, x0) # new initial condition
diff --git a/test/common.jl b/test/common.jl
index 696af3b7e..2c215c2ce 100644
--- a/test/common.jl
+++ b/test/common.jl
@@ -28,19 +28,19 @@ import Logging: Warn
         @test abs(sol[end] - sin(sol.t[end])) < 1e-12
     end
 
-    # @testset "Test integration of ODE with abstract arrays in common interface" begin
-    #     u0 = rand(4, 2)
-    #     tspan = (0.0, 1.0)
-    #     prob = ODEProblem(f!, u0, tspan)
-    #     sol = solve(prob, TaylorMethod(50), abstol=1e-20)
-    #     @test norm(sol[end] - u0.*exp(1)) < 1e-12
-    #     f_oop(u, p, t) = u
-    #     prob_oop = ODEProblem(f_oop, u0, tspan)
-    #     sol_oop = solve(prob_oop, TaylorMethod(50), abstol=1e-20)
-    #     @test norm(sol_oop[end] - u0.*exp(1)) < 1e-12
-    #     @test sol.t == sol_oop.t
-    #     @test sol.u == sol_oop.u
-    # end
+    @testset "Test integration of ODE with abstract arrays in common interface" begin
+        u0 = rand(4, 2)
+        tspan = (0.0, 1.0)
+        prob = ODEProblem(f!, u0, tspan)
+        sol = solve(prob, TaylorMethod(50), abstol=1e-20)
+        @test norm(sol[end] - u0.*exp(1)) < 1e-12
+        @taylorize f_oop(u, p, t) = u
+        prob_oop = ODEProblem(f_oop, u0, tspan)
+        sol_oop = solve(prob_oop, TaylorMethod(50), abstol=1e-20)
+        @test norm(sol_oop[end] - u0.*exp(1)) < 1e-12
+        @test sol.t == sol_oop.t
+        @test sol.u == sol_oop.u
+    end
 
     local tspan = (0.0,5.0)
     local saveat_inputs = ([], 0:1:(tspan[2]+5), 0:1:tspan[2], 3:1:tspan[2], collect(0:1:tspan[2]))
@@ -219,9 +219,14 @@ import Logging: Warn
         @test sol1.alg.parse_eqs == false
         sol2 = (@test_logs min_level=Logging.Warn solve(
             prob, TaylorMethod(order), abstol=abstol))
-        @test sol2.alg.parse_eqs == true
+        sol3 = (@test_logs min_level=Logging.Warn solve(
+            prob, TaylorMethod(order), abstol=abstol, parse_eqs=false))
+        @test sol2.alg.parse_eqs
+        @test !sol3.alg.parse_eqs
         @test sol1.t == sol2.t
         @test sol1.u == sol2.u
+        @test sol1.t == sol3.t
+        @test sol1.u == sol3.u
 
         @taylorize function integ_vec(dx, x, p, t)
             local λ = p[1]
@@ -285,6 +290,81 @@ import Logging: Warn
         @test norm( H_pcr3bp(sol2.u[end]) - J0 ) < 1e-10
         @test sol1.u == sol2.u
         @test sol1.t == sol2.t
+
+        # Same as harmosc1dchain!, but adding a `@threads for`
+        @taylorize function harmosc1dchain!(dq, q, params, t)
+            local N = Int(length(q)/2)
+            local _eltype_q_ = eltype(q)
+            local μ = params
+            X = Array{_eltype_q_}(undef, N, N)
+            accX = Array{_eltype_q_}(undef, N) #acceleration
+            for j in 1:N
+                accX[j] = zero(q[1])
+                dq[j] = q[N+j]
+            end
+            #compute accelerations
+            for j in 1:N
+                for i in 1:N
+                    if i == j
+                    else
+                        X[i,j] = q[i]-q[j]
+                        temp_001 = accX[j] + (μ[i]*X[i,j])
+                        accX[j] = temp_001
+                    end #if i != j
+                end #for, i
+            end #for, j
+            for i in 1:N
+                dq[N+i] = accX[i]
+            end
+            nothing
+        end
+
+        @taylorize function harmosc1dchain_threads!(dq, q, params, t)
+            local N = Int(length(q)/2)
+            local _eltype_q_ = eltype(q)
+            local μ = params
+            X = Array{_eltype_q_}(undef, N, N)
+            accX = Array{_eltype_q_}(undef, N) #acceleration
+            for j in 1:N
+                accX[j] = zero(q[1])
+                dq[j] = q[N+j]
+            end
+            #compute accelerations
+            Threads.@threads for j in 1:N
+                for i in 1:N
+                    if i == j
+                    else
+                        X[i,j] = q[i]-q[j]
+                        temp_001 = accX[j] + (μ[i]*X[i,j])
+                        accX[j] = temp_001
+                    end #if i != j
+                end #for, i
+            end #for, j
+            for i in 1:N
+                dq[N+i] = accX[i]
+            end
+            nothing
+        end
+
+        N = 200
+        x0 = 10randn(2N)
+        μ = 1e-7rand(N)
+
+        @show Threads.nthreads()
+
+        probHO = ODEProblem(harmosc1dchain!, x0, (0.0, 10000.0), μ)
+        solHO = (@test_logs min_level=Logging.Warn solve(
+            probHO, TaylorMethod(order), abstol=abstol, parse_eqs=false))
+        probHOT = ODEProblem(harmosc1dchain_threads!, x0, (0.0, 10000.0), μ)
+        solHOT = (@test_logs min_level=Logging.Warn solve(
+            probHOT, TaylorMethod(order), abstol=abstol))
+        solHOnT = (@test_logs min_level=Logging.Warn solve(
+            probHOT, TaylorMethod(order), abstol=abstol, parse_eqs=false))
+
+        @test solHOT.u == solHOnT.u
+        @test solHOT.t == solHOnT.t
+        @test solHO.u == solHOT.u
+        @test solHO.t == solHOT.t
     end
 
     @testset "Test throwing errors in common interface" begin
diff --git a/test/taylorize.jl b/test/taylorize.jl
index 242d0ec6e..e1b5fcac0 100644
--- a/test/taylorize.jl
+++ b/test/taylorize.jl
@@ -88,24 +88,24 @@ import Logging: Warn
         # TODO: Use metaprogramming here
         tT = t0 + Taylor1(_order)
         xT = x0 + zero(tT)
-        tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
+        rv = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
             Val(xdot1), tT, xT, nothing))
         @test_logs min_level=Logging.Warn TI.jetcoeffs!(
-            Val(xdot1), tT, xT, nothing, tmpTaylor, arrTaylor)
+            Val(xdot1), tT, xT, nothing, rv)
         @test xT ≈ exact_sol(tT, b1, x0)
 
         xT = x0 + zero(tT)
-        tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
+        rv = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
             Val(xdot2), tT, xT, nothing))
         @test_logs min_level=Logging.Warn TI.jetcoeffs!(
-            Val(xdot2), tT, xT, nothing, tmpTaylor, arrTaylor)
+            Val(xdot2), tT, xT, nothing, rv)
         @test xT ≈ exact_sol(tT, 3.0, x0)
 
         xT = x0 + zero(tT)
-        tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
+        rv = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
             Val(xdot2), tT, xT, b1))
         @test_logs min_level=Logging.Warn TI.jetcoeffs!(
-            Val(xdot2), tT, xT, b1, tmpTaylor, arrTaylor)
+            Val(xdot2), tT, xT, b1, rv)
         @test xT ≈ exact_sol(tT, b1, x0)
 
         # The macro returns a (parsed) jetcoeffs! function which yields a `MethodError`
@@ -126,10 +126,10 @@ import Logging: Warn
         @test iszero( norm(tv2-tv2e, Inf) )
         @test iszero( norm(xv2-xv2e, Inf) )
         xT = x0 + zero(tT)
-        tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
+        rv = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
             Val(xdot2_err), tT, xT, nothing))
-        @test_throws MethodError TI.jetcoeffs!(Val(xdot2_err), tT, xT, nothing,
-            tmpTaylor, arrTaylor)
+        @test_throws MethodError TI.jetcoeffs!(
+            Val(xdot2_err), tT, xT, nothing, rv)
 
         # Output includes Taylor polynomial solution
         tv3t, xv3t, polynV3t = (@test_logs (Warn, max_iters_reached()) taylorinteg(
@@ -200,24 +200,24 @@ import Logging: Warn
         # Check that the parsed `jetcoeffs` produces the correct series in `x` and no errors
         tT = 1.0 + Taylor1(_order)
         xT = 10.0 + zero(tT)
-        tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
+        rv = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
             Val(xdot1_parsed), tT, xT, nothing))
         @test_logs min_level=Logging.Warn TI.jetcoeffs!(
-            Val(xdot1_parsed), tT, xT, nothing, tmpTaylor, arrTaylor)
+            Val(xdot1_parsed), tT, xT, nothing, rv)
         @test xT ≈ exact_sol(tT, -10, 10)
 
         xT = 10.0 + zero(tT)
-        tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
+        rv = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
             Val(xdot2), tT, xT, nothing))
         @test_logs min_level=Logging.Warn TI.jetcoeffs!(
-            Val(xdot2), tT, xT, nothing, tmpTaylor, arrTaylor)
+            Val(xdot2), tT, xT, nothing, rv)
         @test xT ≈ exact_sol(tT, -10, 10)
 
         xT = 10.0 + zero(tT)
-        tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
+        rv = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
             Val(xdot3), tT, xT, -10))
         @test_logs min_level=Logging.Warn TI.jetcoeffs!(
-            Val(xdot3), tT, xT, -10, tmpTaylor, arrTaylor)
+            Val(xdot3), tT, xT, -10, rv)
         @test xT ≈ exact_sol(tT, -10, 10)
     end
 
@@ -253,10 +253,10 @@ import Logging: Warn
         tT = t0 + Taylor1(_order)
         qT = q0 .+ zero(tT)
         dqT = similar(qT)
-        tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
+        rv = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
             Val(pendulum!), tT, qT, dqT, nothing))
         @test_logs min_level=Logging.Warn TI.jetcoeffs!(
-            Val(pendulum!), tT, qT, dqT, nothing, tmpTaylor, arrTaylor)
+            Val(pendulum!), tT, qT, dqT, nothing, rv)
     end
 
 
@@ -338,12 +338,12 @@ import Logging: Warn
         tT = t0 + Taylor1(_order)
         zT = zz0 .+ zero(tT)
         dzT = similar(zT)
-        tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
+        rv = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
             Val(eqs3!), tT, zT, dzT, nothing))
         @test_logs min_level=Logging.Warn TI.jetcoeffs!(
-            Val(eqs3!), tT, zT, dzT, nothing, tmpTaylor, arrTaylor)
-        @test typeof(tmpTaylor) == Vector{Taylor1{ComplexF64}}
-        @test typeof(arrTaylor) == Vector{Vector{Taylor1{ComplexF64}}}
+            Val(eqs3!), tT, zT, dzT, nothing, rv)
+        @test typeof(rv.v0) == Vector{Taylor1{ComplexF64}}
+        @test typeof(rv.v1) == Vector{Vector{Taylor1{ComplexF64}}}
         @test zT[1] == exact_sol(tT, -1, z0)
         @test zT[2] == exact_sol(tT, z0, z0)
     end
@@ -404,6 +404,27 @@ import Logging: Warn
         # Compare to exact solution
         @test norm(xv12[end, 1] - sin(tv12[end]), Inf) < 1.0e-15
         @test norm(xv12[end, 2] - cos(tv12[end]), Inf) < 1.0e-15
+
+        @taylorize function fff!(du, u, p, t)
+            x, t = u
+            du[1] = x
+            du[2] = t
+            return du
+        end
+
+        @taylorize function ggg!(du, u, p, t)
+            x, t2 = u
+            du[1] = x
+            du[2] = t2
+            return du
+        end
+
+        tvF, xvF = taylorinteg(ggg!, [1.0, 0.0], 0.0, 10.0, 20, 1.0e-20, parse_eqs=false)
+        tv1, xv1 = taylorinteg(ggg!, [1.0, 0.0], 0.0, 10.0, 20, 1.0e-20)
+        tv2, xv2 = taylorinteg(fff!, [1.0, 0.0], 0.0, 10.0, 20, 1.0e-20)
+
+        @test tv1 == tv2 == tvF
+        @test xv1 == xv2 == xvF
     end
 
 
@@ -438,9 +459,10 @@ import Logging: Warn
         tT = t0 + Taylor1(_order)
         qT = q0 .+ zero(tT)
         dqT = similar(qT)
-        tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
+        rv = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
             Val(harm_osc!), tT, qT, dqT, [1.0]))
-        @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(harm_osc!), tT, qT, dqT, [1.0], tmpTaylor, arrTaylor)
+        @test_logs min_level=Logging.Warn TI.jetcoeffs!(
+            Val(harm_osc!), tT, qT, dqT, [1.0], rv)
         @test qT[1] ≈ cos(tT)
         @test qT[2] ≈ -sin(tT)
 
@@ -466,10 +488,10 @@ import Logging: Warn
         tT = t0 + Taylor1(_order)
         qT = q0 .+ zero(tT)
         dqT = similar(qT)
-        tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
+        rv = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
             Val(harm_osc_error!), tT, qT, dqT, [1.0]))
-        @test_throws MethodError  TI.jetcoeffs!(Val(harm_osc_error!), tT, qT, dqT, [1.0],
-            tmpTaylor, arrTaylor)
+        @test_throws MethodError  TI.jetcoeffs!(
+            Val(harm_osc_error!), tT, qT, dqT, [1.0], rv)
     end
 
     local tf = 100.0
@@ -542,19 +564,22 @@ import Logging: Warn
         tT = 0.0 + Taylor1(_order)
         qT = q0 .+ zero(tT)
         dqT = similar(qT)
-        tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
+        rv = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
             Val(multpendula1!), tT, qT, dqT, (NN, nnrange)))
-        @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(multpendula1!), tT, qT, dqT, (NN, nnrange), tmpTaylor, arrTaylor)
+        @test_logs min_level=Logging.Warn TI.jetcoeffs!(
+            Val(multpendula1!), tT, qT, dqT, (NN, nnrange), rv)
 
         qT = q0 .+ zero(tT)
-        tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
+        rv = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
             Val(multpendula2!), tT, qT, dqT, nothing))
-        @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(multpendula2!), tT, qT, dqT, nothing, tmpTaylor, arrTaylor)
+        @test_logs min_level=Logging.Warn TI.jetcoeffs!(
+            Val(multpendula2!), tT, qT, dqT, nothing, rv)
 
         qT = q0 .+ zero(tT)
-        tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
+        rv = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
             Val(multpendula3!), tT, qT, dqT, [NN, nnrange]))
-        @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(multpendula3!), tT, qT, dqT, [NN, nnrange], tmpTaylor, arrTaylor)
+        @test_logs min_level=Logging.Warn TI.jetcoeffs!(
+            Val(multpendula3!), tT, qT, dqT, [NN, nnrange], rv)
     end
 
 
@@ -671,27 +696,31 @@ import Logging: Warn
         tT = 0.0 + Taylor1(_order)
         qT = q0 .+ zero(tT)
         dqT = similar(qT)
-        tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(Val(kepler1!),
-            tT, qT, dqT, -1.0))
-        @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(kepler1!), tT, qT, dqT, -1.0, tmpTaylor, arrTaylor)
+        rv = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
+            Val(kepler1!), tT, qT, dqT, -1.0))
+        @test_logs min_level=Logging.Warn TI.jetcoeffs!(
+            Val(kepler1!), tT, qT, dqT, -1.0, rv)
 
         qT = q0 .+ zero(tT)
         dqT = similar(qT)
-        tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(Val(kepler2!),
-            tT, qT, dqT, pars))
-        @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(kepler2!), tT, qT, dqT, pars, tmpTaylor, arrTaylor)
+        rv = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
+            Val(kepler2!), tT, qT, dqT, pars))
+        @test_logs min_level=Logging.Warn TI.jetcoeffs!(
+            Val(kepler2!), tT, qT, dqT, pars, rv)
 
         qT = q0 .+ zero(tT)
         dqT = similar(qT)
-        tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(Val(kepler3!),
-            tT, qT, dqT, nothing))
-        @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(kepler3!), tT, qT, dqT, nothing, tmpTaylor, arrTaylor)
+        rv = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
+            Val(kepler3!), tT, qT, dqT, nothing))
+        @test_logs min_level=Logging.Warn TI.jetcoeffs!(
+            Val(kepler3!), tT, qT, dqT, nothing, rv)
 
         qT = q0 .+ zero(tT)
         dqT = similar(qT)
-        tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(Val(kepler4!),
-            tT, qT, dqT, nothing))
-        @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(kepler4!), tT, qT, dqT, nothing, tmpTaylor, arrTaylor)
+        rv = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
+            Val(kepler4!), tT, qT, dqT, nothing))
+        @test_logs min_level=Logging.Warn TI.jetcoeffs!(
+            Val(kepler4!), tT, qT, dqT, nothing, rv)
     end
 
 
@@ -809,24 +838,28 @@ import Logging: Warn
         tT = 0.0 + Taylor1(_order)
         qT = q0 .+ zero(tT)
         dqT = similar(qT)
-        tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(Val(kepler1!),
-            tT, qT, dqT, nothing))
-        @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(kepler1!), tT, qT, dqT, nothing, tmpTaylor, arrTaylor)
+        rv = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
+            Val(kepler1!), tT, qT, dqT, nothing))
+        @test_logs min_level=Logging.Warn TI.jetcoeffs!(
+            Val(kepler1!), tT, qT, dqT, nothing, rv)
 
         qT = q0 .+ zero(tT)
-        tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(Val(kepler2!),
-            tT, qT, dqT, -1.0))
-        @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(kepler2!), tT, qT, dqT, -1.0, tmpTaylor, arrTaylor)
+        rv = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
+            Val(kepler2!), tT, qT, dqT, -1.0))
+        @test_logs min_level=Logging.Warn TI.jetcoeffs!(
+            Val(kepler2!), tT, qT, dqT, -1.0, rv)
 
         qT = q0 .+ zero(tT)
-        tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(Val(kepler3!),
-            tT, qT, dqT, -1.0))
-        @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(kepler3!), tT, qT, dqT, -1.0, tmpTaylor, arrTaylor)
+        rv = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
+            Val(kepler3!), tT, qT, dqT, -1.0))
+        @test_logs min_level=Logging.Warn TI.jetcoeffs!(
+            Val(kepler3!), tT, qT, dqT, -1.0, rv)
 
         qT = q0 .+ zero(tT)
-        tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(Val(kepler4!),
-            tT, qT, dqT, nothing))
-        @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(kepler4!), tT, qT, dqT, nothing, tmpTaylor, arrTaylor)
+        rv = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
+            Val(kepler4!), tT, qT, dqT, nothing))
+        @test_logs min_level=Logging.Warn TI.jetcoeffs!(
+            Val(kepler4!), tT, qT, dqT, nothing, rv)
     end
 
 
@@ -977,17 +1010,17 @@ import Logging: Warn
         tT = 0.0 + Taylor1(_order)
         qT = q0 .+ zero(tT)
         dqT = similar(qT)
-        tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
+        rv = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
             Val(lorenz1!), tT, qT, dqT, params))
-        @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(lorenz1!),
-            tT, qT, dqT, params, tmpTaylor, arrTaylor)
+        @test_logs min_level=Logging.Warn TI.jetcoeffs!(
+            Val(lorenz1!), tT, qT, dqT, params, rv)
 
         qT = q0 .+ zero(tT)
         dqT = similar(qT)
-        tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
+        rv = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(
             Val(lorenz2!), tT, qT, dqT, nothing))
-        @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(lorenz2!),
-            tT, qT, dqT, nothing, tmpTaylor, arrTaylor)
+        @test_logs min_level=Logging.Warn TI.jetcoeffs!(
+            Val(lorenz2!), tT, qT, dqT, nothing, rv)
     end
 
 
@@ -1032,7 +1065,7 @@ import Logging: Warn
         @test_throws ArgumentError TI._make_parsed_jetcoeffs(ex)
 
         # The macro works fine, but the `jetcoeffs!` method does not work
-        # (so the integration would silently run using `parse_eqs=false`)
+        # (so the integration would run using `parse_eqs=false`)
         @taylorize function harm_osc!(dx, x, p, t)
             local ω = p[1]
             # local ω2 = ω^2    # Needed this to avoid an error
@@ -1042,8 +1075,10 @@ import Logging: Warn
         end
         tT = t0 + Taylor1(_order)
         qT = [1.0, 0.0] .+ zero(tT)
-        @test_throws MethodError TI.__jetcoeffs!(
-            Val(true), harm_osc!, tT, qT, similar(qT), similar(qT), [2.0])
+        parse_eqs, rv = (@test_logs (Warn, unable_to_parse(harm_osc!)) TI._determine_parsing!(
+            true, harm_osc!, tT, qT, similar(qT), [2.0]))
+        @test !parse_eqs
+        @test_throws MethodError TI.__jetcoeffs!(Val(harm_osc!), tT, qT, similar(qT), [2.0], rv)
 
         @taylorize function kepler1!(dq, q, p, t)
             μ = p
@@ -1058,8 +1093,25 @@ import Logging: Warn
         end
         tT = t0 + Taylor1(_order)
         qT = [0.2, 0.0, 0.0, 3.0] .+ zero(tT)
-        @test_throws MethodError TI.__jetcoeffs!(
-            Val(true), kepler1!, tT, qT, similar(qT), similar(qT), -1.0)
+        parse_eqs, rv = (@test_logs (Warn, unable_to_parse(kepler1!)) TI._determine_parsing!(
+            true, kepler1!, tT, qT, similar(qT), -1.0))
+        @test !parse_eqs
+        @test_throws MethodError TI.__jetcoeffs!(Val(kepler1!), tT, qT, similar(qT), -1.0, rv)
+
+        # Error: `@taylorize` allows only to parse up tp 5-index arrays
+        ex = :(function err_arr_indx!(Dz, z, p, t)
+            n = size(z,1)
+            arr3 = Array{typeof(z[1])}(undef, n, 1, 1)
+            arr4 = Array{typeof(z[1])}(undef, n, 1, 1, 1)
+            arr5 = Array{typeof(z[1])}(undef, n, 1, 1, 1, 1)
+            for i in eachindex(z)
+                arr5[i,1,1,1,1] = zero(z[1])
+                Dz[i] = z[i] + arr5[i,1,1,1,1]
+            end
+            nothing
+        end)
+        @test_throws ErrorException("Error: `@taylorize` allows only to parse up tp 5-index arrays") TI._make_parsed_jetcoeffs(ex)
+
     end
 
 
@@ -1207,21 +1259,28 @@ import Logging: Warn
         @test newex1.args[1] == :(
             TaylorIntegration.jetcoeffs!(::Val{f!}, t::Taylor1{_T},
                 q::AbstractArray{Taylor1{_S}, _N},
-                dq::AbstractArray{Taylor1{_S}, _N}, p, __tmpTaylor, __arrTaylor) where
+                dq::AbstractArray{Taylor1{_S}, _N}, p,
+                __ralloc::TaylorIntegration.RetAlloc{Taylor1{_S}}) where
                     {_T <: Real, _S <: Number, _N})
 
         # Include not recognized functions as they appear
-        @test newex1.args[2].args[2] == :(__tmpTaylor[1] = my_simple_function(q, p, t))
+        @test newex1.args[2].args[2] == :(aa = __ralloc.v0[1])
+        @test newex1.args[2].args[3] == :(aa = my_simple_function(q, p, t))
         @test newex2.args[2].args[2] == :(aa = my_simple_function(q, p, t))
-        @test newex1.args[2].args[5].args[2].args[2] ==
-            :(__tmpTaylor[1] = my_simple_function(q, p, t))
+        @test newex1.args[2].args[6].args[2].args[2] ==
+            :(aa = my_simple_function(q, p, t))
 
         # Return line
         @test newex1.args[2].args[end] == :(return nothing)
-        @test newex2.args[2].args[end] == :(return ([aa], [Vector{Taylor1{_S}}(undef, 0)]))
+        @test newex2.args[2].args[end] ==
+            :(return TaylorIntegration.RetAlloc{Taylor1{_S}}([aa],
+                [Array{Taylor1{_S},1}(undef, 0)],
+                [Array{Taylor1{_S},2}(undef, 0, 0)],
+                [Array{Taylor1{_S},3}(undef, 0, 0, 0)],
+                [Array{Taylor1{_S},4}(undef, 0, 0, 0, 0)]))
 
         # Issue 96: deal with `elseif`s, `continue` and `break`
-        @test newex1.args[2].args[5].args[2].args[3] == :(
+        @test newex1.args[2].args[6].args[2].args[3] == :(
             for i = 1:length(q)
                 if i == 1
                     TaylorSeries.mul!(dq[i], 2, q[i], ord)
@@ -1230,7 +1289,7 @@ import Logging: Warn
                         TaylorSeries.identity!(dq[i], q[i], ord)
                     else
                         if i == 3
-                            TaylorSeries.identity!(dq[i], __tmpTaylor[1], ord)
+                            TaylorSeries.identity!(dq[i], aa, ord)
                             continue
                         else
                             dq[i] = my_complicate_function(q)
@@ -1239,6 +1298,27 @@ import Logging: Warn
                     end
                 end
             end)
+
+        # Throws no error
+        ex = :(
+            function err_arr_indx!(Dz, z, p, t)
+                local n = size(z,1)  # important to include it!
+                arr3 = Array{typeof(z[1])}(undef, n, 1, 1)
+                arr4 = Array{typeof(z[1])}(undef, n, 1, 1, 1)
+                for i in eachindex(z)
+                    arr3[i,1,1] = zero(z[1])
+                    arr4[i,1,1,1] = zero(z[1])
+                    Dz[i] = z[i] + arr4[i,1,1,1]
+                end
+                nothing
+            end)
+        newex1, newex2 = TI._make_parsed_jetcoeffs(ex)
+        @test newex1.args[2].args[2] == :(arr3 = __ralloc.v3[1])
+        @test newex1.args[2].args[3] == :(arr4 = __ralloc.v4[1])
+        @test newex2.args[2].args[end].args[1].args[3] == :([Array{Taylor1{_S}, 1}(undef, 0)])
+        @test newex2.args[2].args[end].args[1].args[4] == :([Array{Taylor1{_S}, 2}(undef, 0, 0)])
+        @test newex2.args[2].args[end].args[1].args[5] == :([arr3])
+        @test newex2.args[2].args[end].args[1].args[6] == :([arr4])
     end
 
     @testset "Test @taylorize with @threads" begin
@@ -1288,26 +1368,26 @@ import Logging: Warn
         t3 = deepcopy(t1)
 
         # No error is thrown
-        parse_eqs, tmpTaylor, arrTaylor = (@test_logs min_level=Warn TI._determine_parsing!(
+        parse_eqs, rv = (@test_logs min_level=Warn TI._determine_parsing!(
             true, f1!, t1, x01, dx01, nothing))
         @test parse_eqs
-        @test typeof(tmpTaylor) == Vector{Taylor1{Float64}}
-        @test typeof(arrTaylor) == Vector{Vector{Taylor1{Float64}}}
-        parse_eqs, tmpTaylor, arrTaylor = (@test_logs min_level=Warn TI._determine_parsing!(
+        @test typeof(rv.v0) == Vector{Taylor1{Float64}}
+        @test typeof(rv.v1) == Vector{Vector{Taylor1{Float64}}}
+        parse_eqs, rv = (@test_logs min_level=Warn TI._determine_parsing!(
             true, f1_parsed!, t1p, x01p, dx01p, nothing))
         @test parse_eqs
-        @test typeof(tmpTaylor) == Vector{Taylor1{Float64}}
-        @test typeof(arrTaylor) == Vector{Vector{Taylor1{Float64}}}
-        parse_eqs, tmpTaylor, arrTaylor = (@test_logs min_level=Warn TI._determine_parsing!(
+        @test typeof(rv.v0) == Vector{Taylor1{Float64}}
+        @test typeof(rv.v1) == Vector{Vector{Taylor1{Float64}}}
+        parse_eqs, rv = (@test_logs min_level=Warn TI._determine_parsing!(
             true, f2!, t2, x02, dx02, nothing))
         @test parse_eqs
-        @test typeof(tmpTaylor) == Vector{Taylor1{Float64}}
-        @test typeof(arrTaylor) == Vector{Vector{Taylor1{Float64}}}
-        parse_eqs, tmpTaylor, arrTaylor = (@test_logs min_level=Warn TI._determine_parsing!(
+        @test typeof(rv.v0) == Vector{Taylor1{Float64}}
+        @test typeof(rv.v1) == Vector{Vector{Taylor1{Float64}}}
+        parse_eqs, rv = (@test_logs min_level=Warn TI._determine_parsing!(
             true, f3!, t3, x03, dx03, nothing))
         @test parse_eqs
-        @test typeof(tmpTaylor) == Vector{Taylor1{Float64}}
-        @test typeof(arrTaylor) == Vector{Vector{Taylor1{Float64}}}
+        @test typeof(rv.v0) == Vector{Taylor1{Float64}}
+        @test typeof(rv.v1) == Vector{Vector{Taylor1{Float64}}}
 
         @test x01 == x01p
         @test x01p == x02
@@ -1381,18 +1461,20 @@ import Logging: Warn
 
         @show Threads.nthreads()
 
-        parse_eqs, _, _ = (@test_logs min_level=Warn TI._determine_parsing!(
+        parse_eqs, rv = (@test_logs min_level=Warn TI._determine_parsing!(
             true, harmosc1dchain!, t, x, dx, μ))
         @test parse_eqs
-        parse_eqs_, _, _ = (@test_logs min_level=Warn TI._determine_parsing!(
-            true, harmosc1dchain!, t, x_, dx_, μ))
+        parse_eqs, rv = (@test_logs min_level=Warn TI._determine_parsing!(
+            true, harmosc1dchain_threads!, t, x_, dx_, μ))
         @test parse_eqs
 
         @test x == x_
         @test dx == dx_
 
-        tv, xv = (@test_logs min_level=Warn taylorinteg(harmosc1dchain!, x0, t0, 10000.0, _order, _abstol, μ))
-        tv_, xv_ = (@test_logs min_level=Warn taylorinteg(harmosc1dchain_threads!, x0, t0, 10000.0, _order, _abstol, μ))
+        tv, xv = (@test_logs min_level=Warn taylorinteg(
+            harmosc1dchain!, x0, t0, 10000.0, _order, _abstol, μ))
+        tv_, xv_ = (@test_logs min_level=Warn taylorinteg(
+            harmosc1dchain_threads!, x0, t0, 10000.0, _order, _abstol, μ))
 
         @test tv == tv_
         @test xv == xv_