Update docs

docs/src/api.md

@@ -14,6 +14,12 @@ lyap_taylorinteg
 @taylorize
 ```
 
+## Exported types
+
+```@docs
+TaylorSolution
+```
+
 ## Internal
 
 ```@autodocs

docs/src/lorenz_lyapunov.md

@@ -114,12 +114,14 @@ one with the values of the Lyapunov spectrum.
 
 We first carry out the integration computing internally the Jacobian
 ```@example lorenz
-tv, xv, λv = lyap_taylorinteg(lorenz!, x0, t0, tmax, 28, 1e-20, params; maxsteps=2000000);
+sol = lyap_taylorinteg(lorenz!, x0, t0, tmax, 28, 1e-20, params; maxsteps=2000000);
+tv, xv, λv = sol.t, sol.x, sol.λ
 nothing # hide
 ```
 Now, the integration is obtained exploiting `lorenz_jac!`:
 ```@example lorenz
-tv_, xv_, λv_ = lyap_taylorinteg(lorenz!, x0, t0, tmax, 28, 1e-20, params, lorenz_jac!; maxsteps=2000000);
+sol_ = lyap_taylorinteg(lorenz!, x0, t0, tmax, 28, 1e-20, params, lorenz_jac!; maxsteps=2000000);
+tv_, xv_, λv_ = sol_.t, sol_.x, sol_.λ
 nothing # hide
 ```
 In terms of performance the second method is about ~50% faster than the first.

docs/src/pendulum.md

@@ -103,7 +103,7 @@ We perform the Taylor integration using the initial condition `x0TN`, during
 6 periods of the pendulum (i.e., ``6T``), exploiting multiple dispatch:
 ```@example pendulum
 tv = t0:integstep:tmax # the times at which the solution will be evaluated
-xv = taylorinteg(pendulum!, q0TN, tv, order, abstol)
+sol = taylorinteg(pendulum!, q0TN, tv, order, abstol)
 nothing # hide
 ```
 
@@ -128,19 +128,19 @@ nothing # hide
 ```
 
 We evaluate the jet at ``\partial U_{x(t)}`` (the boundary of ``U_{x(t)}``) at each
-value of the solution vector `xv`; we organize these values such that we can plot
+value of the solution array `sol.x`; we organize these values such that we can plot
 them later:
 ```@example pendulum
-xjet_plot = map(λ->λ.(ξv), xv[:,1])
-pjet_plot = map(λ->λ.(ξv), xv[:,2])
+xjet_plot = map(λ->λ.(ξv), sol.x[:,1])
+pjet_plot = map(λ->λ.(ξv), sol.x[:,2])
 nothing # hide
 ```
 Above, we have exploited the fact that `Array{TaylorN{Float64}}` variables are
 callable objects. Now, we evaluate the jet at the nominal solution, which
-corresponds to ``\xi=(0,0)``, at each value of the solution vector `xv`:
+corresponds to ``\xi=(0,0)``, at each value of the solution vector `sol.x`:
 ```@example pendulum
-x_nom = xv[:,1]()
-p_nom = xv[:,2]()
+x_nom = sol.x[:,1]()
+p_nom = sol.x[:,2]()
 nothing # hide
 ```
 

docs/src/root_finding.md

@@ -124,11 +124,11 @@ for i in 1:nconds
     x_ini .= x0 .+ 0.005 .* [0.0, sqrt(rand1)*cos(2pi*rand2), 0.0, sqrt(rand1)*sin(2pi*rand2)]
     px!(x_ini, E0)   # ensure initial energy is E0
 
-    tv_i, xv_i, tvS_i, xvS_i, gvS_i = taylorinteg(henonheiles!, g, x_ini, 0.0, 135.0,
+    sol_i = taylorinteg(henonheiles!, g, x_ini, 0.0, 135.0,
         25, 1e-25, maxsteps=30000);
-    tvSv[i] = vcat(0.0, tvS_i)
-    xvSv[i] = vcat(transpose(x_ini), xvS_i)
-    gvSv[i] = vcat(0.0, gvS_i)
+    tvSv[i] = vcat(0.0, sol_i.t)
+    xvSv[i] = vcat(transpose(x_ini), sol_i.x)
+    gvSv[i] = vcat(0.0, sol_i.gresids)
 end
 nothing # hide
 ```
@@ -202,14 +202,15 @@ nothing # hide
 
 We are now set to carry out the integration.
 ```@example poincare
-tvTN, xvTN, tvSTN, xvSTN, gvSTN = taylorinteg(henonheiles!, g, x0TN, 0.0, 135.0, 25, 1e-25, maxsteps=30000);
+solTN = taylorinteg(henonheiles!, g, x0TN, 0.0, 135.0, 25, 1e-25, maxsteps=30000);
 nothing # hide
 ```
 
 We define some auxiliary arrays, and then make an animation with the results for
 plotting.
 ```@example poincare
 #some auxiliaries:
+tvTN, xvTN, tvSTN, xvSTN, gvSTN = solTN.t, solTN.x, solTN.tevents, solTN.xevents, solTN.gresids
 xvSTNaa = Array{Array{TaylorN{Float64},1}}(undef, length(tvSTN)+1 );
 xvSTNaa[1] = x0TN
 for ind in 2:length(tvSTN)+1

docs/src/simple_example.md

@@ -112,7 +112,7 @@ but it actually does not reach this value.
 Figure 1 displays the computed solution as a function of
 time, in log scale.
 ```@example example1
-plot(sol.t, log10.sol.x), shape=:circle)
+plot(sol.t, log10.(sol.x), shape=:circle)
 xlabel!("t")
 ylabel!("log10(x(t))")
 xlims!(0,0.34)

docs/src/taylorize.md

@@ -56,7 +56,7 @@ tf = 10000.0
 q0 = [1.3, 0.0]
 
 # The actual integration
-t1, x1 = taylorinteg(pendulumNP!, q0, t0, tf, 25, 1e-20, maxsteps=50000); # warm-up run
+sol1 = 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
@@ -119,7 +119,7 @@ use the same instruction as above; the default value for the keyword argument `p
 is `true`, so we may omit it.
 
 ```@example taylorize
-t2, x2 = taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=50000); # warm-up run
+sol2 = 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
@@ -132,7 +132,7 @@ e1/e2, all1/all2
 
 We can check that both integrations yield the same results.
 ```@example taylorize
-t1 == t2 && x1 == x2
+sol1.t == sol2.t && sol1.x == sol2.x
 ```
 
 As stated above, in order to allow opting out of using the specialized method

src/integrator/taylorsolution.jl

@@ -4,6 +4,15 @@ abstract type AbstractTaylorSolution{T<:Real, U<:Number} end
 
 ## Constructors
 
+"""
+TaylorSolution{T, U, N, VT<:AbstractVector{T}, AX<:AbstractArray{U,N},
+        P<:Union{Nothing, AbstractArray{Taylor1{U}, N}},
+        VTE<:Union{Nothing, AbstractVector{U}},
+        AXE<:Union{Nothing, AbstractArray{U, N}},
+        VΛ<:Union{Nothing, AbstractArray{U,N}}} <: AbstractTaylorSolution{T, U}
+
+Return type for `taylorinteg`.
+"""
 struct TaylorSolution{T, U, N, VT<:AbstractVector{T}, AX<:AbstractArray{U,N},
         P<:Union{Nothing, AbstractArray{Taylor1{U}, N}},
         VTE<:Union{Nothing, AbstractVector{U}},