diff --git a/src/vfvm_solver.jl b/src/vfvm_solver.jl
index 6120430e3..e96c6394f 100644
--- a/src/vfvm_solver.jl
+++ b/src/vfvm_solver.jl
@@ -397,13 +397,13 @@ function CommonSolve.solve(
                 end
                 if solved
                     Δu = control.delta(system, solution, oldsolution, λ, Δλ)
-                    if Δu > 2.0 * control.Δu_opt
+                    if Δu > control.Δu_max_factor * control.Δu_opt
                         solved = false
                     end
                 end
                 if !solved
                     # reduce time step and retry  solution
-                    Δλ = Δλ * 0.5
+                    Δλ = Δλ * control.Δt_decrease
                     rd(x)=round(x,sigdigits=5)
                     if Δλ < Δλ_min(control, transient)
                         if !(control.force_first_step && istep == 0)
@@ -478,6 +478,7 @@ function CommonSolve.solve(
     end
 
     system.history = allhistory
+    tsol.history = allhistory
     return tsol
 end
 
@@ -486,8 +487,8 @@ end
                                    t=0.0, tstep=Inf,embed=0.0)
 
 Evaluate residual and jacobian at solution value u.
-Returns a solution vector containing the residual, and an ExendableSparseMatrix
-containing the linearization at u.
+Returns a solution vector containing a copy of  residual, and an ExendableSparseMatrix
+containing a copy of the linearization at u.
 
 """
 function evaluate_residual_and_jacobian(sys, u; t = 0.0, tstep = Inf, embed = 0.0)
diff --git a/src/vfvm_solvercontrol.jl b/src/vfvm_solvercontrol.jl
index e9b951a95..3a34db181 100644
--- a/src/vfvm_solvercontrol.jl
+++ b/src/vfvm_solvercontrol.jl
@@ -184,12 +184,23 @@ Base.@kwdef mutable struct SolverControl
     Δt_grow::Float64 = 1.2
 
     """
-    Optimal size of update for time stepping and embeding.
+    Time step decrease factor upon rejection
+    """
+    Δt_decrease::Float64 = 0.5
+
+    """
+    Optimal size of update for time stepping and embedding.
     The algorithm tries to keep the difference in norm between "old" and "new" 
     solutions  approximately at this value.
     """
     Δu_opt::Float64 = 0.1
 
+    """
+    Control maximum  sice of update `Δu` for time stepping and embeding relative to
+    `Δu_opt`. Time steps with `Δu > Δu_max_factor*Δu_opt` will be rejected.
+    """
+    Δu_max_factor::Float64 = 1.2
+
     """
     Force first timestep.
     """
@@ -242,7 +253,8 @@ Base.@kwdef mutable struct SolverControl
     sample::Function = function (sol, t) end
 
     """
-    Time step error estimator
+    Time step error estimator. A function `Δu=delta(system,u,uold,t,Δt)` to
+    calculate `Δu`.
     """
     delta::Function = (system, u, v, t, Δt) -> norm(system, u - v, Inf)
 
diff --git a/src/vfvm_transientsolution.jl b/src/vfvm_transientsolution.jl
index a545f72b9..5d8a010ac 100644
--- a/src/vfvm_transientsolution.jl
+++ b/src/vfvm_transientsolution.jl
@@ -33,22 +33,54 @@ mutable struct TransientSolution{T, N, A, B} <: AbstractTransientSolution{T, N,
     Vector of solutions
     """
     u::A
+
     """
     Vector of times
     """
     t::B
+
+    """
+    History
+    """
+    history::TransientSolverHistory
 end
 
 function TransientSolution(vec::AbstractVector{T}, ts, ::NTuple{N}) where {T, N}
-    TransientSolution{T, N, typeof(vec), typeof(ts)}(vec, ts)
+    TransientSolution{T, N, typeof(vec), typeof(ts)}(vec, ts, TransientSolverHistory())
 end
 
+
 TransientSolution(vec::AbstractVector, ts::AbstractVector) = TransientSolution(vec, ts, (size(vec[1])..., length(vec)))
 
 Base.append!(s::AbstractTransientSolution, t::Real, sol::AbstractArray) = push!(s.t, t), push!(s.u, copy(sol))
 
 (sol::AbstractTransientSolution)(t) = _interpolate(sol, t)
 
+
+"""
+    history(tsol)
+
+Return solver history if `log` was set to true.
+See  see [`NewtonSolverHistory`](@ref), [`TransientSolverHistory`](@ref).
+"""
+history(tsol::AbstractTransientSolution) = tsol.history
+
+"""
+    history_details(tsol)
+
+Return details of solver history from last `solve` call, if `log` was set to true.
+See [`details`](@ref).
+"""
+history_details(tsol::AbstractTransientSolution) = details(sys.history)
+
+"""
+    history_summary(tsol)
+
+Return summary of solver history from last `solve` call, if `log` was set to true.
+"""
+history_summary(tsol::AbstractTransientSolution) = summary(tsol.history)
+
+
 function _interpolate(sol, t)
     if isapprox(t, sol.t[1]; atol = 1.0e-10 * abs(sol.t[2] - sol.t[1]))
         return sol[1]