diff --git a/lib/ControlSystemsBase/src/timeresp.jl b/lib/ControlSystemsBase/src/timeresp.jl
index 806e9a8ac..e28e871cc 100644
--- a/lib/ControlSystemsBase/src/timeresp.jl
+++ b/lib/ControlSystemsBase/src/timeresp.jl
@@ -138,7 +138,7 @@ plot(result, plotu=true, plotx=false)
 ```
 `y`, `x`, `u` have time in the second dimension. Initial state `x0` defaults to zero.
 
-Continuous-time systems are simulated using an ODE solver if `u` is a function (requires using ControlSystems). If `u` is an array, the system is discretized (with `method=:zoh` by default) before simulation. For a lower-level interface, see `?Simulator` and `?solve`. For continuous-time systems, keyword arguments are forwarded to the ODE solver. By default, the option `dtmax = t[2]-t[1]` is used to prevent the solver from stepping over discontinuities in `u(x, t)`. This prevents the solver from taking too large steps, but may also slow down the simulation when `u` is smooth. To disable this behavior, set `dtmax = Inf`.
+Continuous-time systems are simulated using an ODE solver if `u` is a function (requires using ControlSystems). If `u` is an array, the system is discretized (with `method=:zoh` by default) before simulation. For a lower-level interface, see `?Simulator` and `?solve`. For continuous-time systems, keyword arguments are forwarded to the ODE solver. By default, the option `dtmax = t[2]-t[1]` is used to prevent the solver from stepping over discontinuities in `u(x, t)`. This prevents the solver from taking too large steps, but may also slow down the simulation when `u` is smooth. To disable this behavior, set `dtmax = Inf`. The default ODE solver is `AutoTsit5(Rosenbrock23())`, which switches from the explicit method `Tsit5` to the method `Rosenbrock23` for stiff problems when stiffness is detected. `Rosenbrock23` computes derivatives by forward-mode automatic differentiation and therefore evaluates `u` with arguments of a dual-number type; if `u` does not support such arguments, select another solver with the keyword argument `alg`, e.g., `alg = Tsit5()`.
 
 `u` can be a function or a *matrix* of precalculated control signals and must have dimensions `(nu, length(t))`.
 If `u` is a function, then `u(x,i)` (for discrete systems) or `u(x,t)` (for continuous ones) is called to calculate the control signal at every iteration (time instance used by solver). This can be used to provide a control law such as state feedback `u(x,t) = -L*x` calculated by `lqr`.
diff --git a/src/timeresp.jl b/src/timeresp.jl
index 75156de8d..c11c041d7 100644
--- a/src/timeresp.jl
+++ b/src/timeresp.jl
@@ -1,4 +1,4 @@
-import OrdinaryDiffEq: ODEProblem, Tsit5, solve
+import OrdinaryDiffEq: ODEProblem, Tsit5, AutoTsit5, Rosenbrock23, solve
 import OrdinaryDiffEqLowOrderRK: BS3
 import ControlSystemsBase: lsim, step, impulse, HammersteinWienerSystem, DelayLtiSystem, SimResult
 import DelayDiffEq: MethodOfSteps
@@ -21,9 +21,19 @@ Internal function: Dynamics equation for simulation of a linear system.
     mul!(dx, B, u(x, t), 1, 1)
 end
 
+"""
+    check_retcode(sol)
+
+Internal function: Throw an error if the differential-equation solver did not complete the simulation successfully.
+"""
+function check_retcode(sol)
+    OrdinaryDiffEq.SciMLBase.successful_retcode(sol) && return
+    error("The simulation terminated at t = $(sol.t[end]) with return code $(sol.retcode). Consider selecting another solver with the keyword argument `alg`, e.g., a solver for stiff problems, or adjusting the keyword arguments `abstol`, `reltol` and `maxiters`.")
+end
+
 # This method is more specific than the lsim in ControlSystemsBase that does not specify Continuous timeevol for sys, hence, if ControlSystems is loaded, ControlSystems.lsim will take precedence over ControlSystemsBase.lsim for Continuous systems
 function lsim(sys::AbstractStateSpace{Continuous}, u::Function, t::AbstractVector;
-        x0::AbstractVecOrMat=zeros(Bool, nstates(sys)), method::Symbol=:cont, dtmax = t[2]-t[1], alg = Tsit5(), kwargs...)
+        x0::AbstractVecOrMat=zeros(Bool, nstates(sys)), method::Symbol=:cont, dtmax = t[2]-t[1], alg = AutoTsit5(Rosenbrock23()), kwargs...)
     ny, nu = size(sys)
     nx = sys.nx
     u0 = u(x0,t[1])
@@ -56,6 +66,7 @@ function lsim(sys::AbstractStateSpace{Continuous}, u::Function, t::AbstractVecto
     else
         p = (sys.A, sys.B, u)
         sol = solve(ODEProblem(f_lsim, x0, (t[1], t[end]+dt/2), p), alg; dtmax, saveat=t, kwargs...)
+        check_retcode(sol)
         x = reduce(hcat, sol.u)::Matrix{T}
         uout = Matrix{T}(undef, nu, length(t))
         for i = eachindex(t)
@@ -216,6 +227,7 @@ function _lsim(sys::DelayLtiSystem{T,S}, Base.@nospecialize(u!), t::AbstractArra
                 callback=cb)
     # Important to stop at t since we can not access derivatives in SavingCallback
     sol = DelayDiffEq.solve(prob, alg; tstops=t, saveat=t, kwargs...)
+    check_retcode(sol)
 
     # Retrive the saved values
     uout2 = zeros(T, nu, nt)
@@ -412,6 +424,7 @@ function _lsim(sys::HammersteinWienerSystem{T}, u!, t::AbstractArray{<:Real}, x0
         p = (A, B1, B2, C1, C2, D11, D12, D21, D22, f, u!, uout, dy, order)
         prob = ODEProblem{true}(hw_f, x0, (T(t[1]), T(t[end]+dt/2)), p)
         sol = OrdinaryDiffEq.solve(prob, alg; saveat=t, kwargs...)
+        check_retcode(sol)
         x = reduce(hcat, sol.u)::Matrix{T}
     else
         x = zeros(T, 0, nt)        # Empty State matrix
diff --git a/test/test_delay_timeresp.jl b/test/test_delay_timeresp.jl
index 8e010719c..afa5f4e85 100644
--- a/test/test_delay_timeresp.jl
+++ b/test/test_delay_timeresp.jl
@@ -144,3 +144,6 @@ let sys = delay(0.01)
     @test res.y[1] ≈ 0 atol = 1e-8
     @test res.y[end] ≈ 1 atol = 1e-3
 end
+
+# An unsuccessful simulation throws an error that states the return code
+@test_throws "return code" lsim(delay(1)*tf(1, [1, 1]), (x, t) -> [1.0], 0:0.1:10; maxiters=5)
diff --git a/test/test_nonlinear_timeresp.jl b/test/test_nonlinear_timeresp.jl
index 633998693..b0230b3dc 100644
--- a/test/test_nonlinear_timeresp.jl
+++ b/test/test_nonlinear_timeresp.jl
@@ -257,3 +257,6 @@ res_smooth = lsim(sys_hyst_smooth, (x,t) -> 5sin(t), t)
 @test maximum(res_smooth.y) ≈ amplitude rtol=0.15
 @test minimum(res_smooth.y) ≈ -amplitude rtol=0.15
 
+
+# An unsuccessful simulation throws an error that states the return code
+@test_throws "return code" lsim(feedback(saturation(0.5)*tf(1, [1, 1])), (x, t) -> [1.0], 0:0.1:10; maxiters=5)
diff --git a/test/test_timeresp.jl b/test/test_timeresp.jl
index c4fba6f44..b59062f83 100644
--- a/test/test_timeresp.jl
+++ b/test/test_timeresp.jl
@@ -49,6 +49,17 @@ uD(x, t) = [3.0]
 resD = lsim(sysD, uD, tD)
 @test resD.y[:] == fill(2.0 * 3.0, length(tD))
 
+# Stiff system with one slow and one fast pole
+sysstiff = ss(tf(1, [1, 1e6+1, 1e6]))
+tstiff = 0:0.1:10
+ustiff(x, t) = [1.0]
+resstiff = lsim(sysstiff, ustiff, tstiff)
+reszoh = lsim(sysstiff, ones(1, length(tstiff)), tstiff) # Exact for a constant input
+@test size(resstiff.y) == (1, length(tstiff))
+@test resstiff.y ≈ reszoh.y rtol=1e-3
+# An unsuccessful simulation throws an error that states the return code
+@test_throws "return code" lsim(sysstiff, ustiff, tstiff; alg=OrdinaryDiffEq.Tsit5(), maxiters=100)
+
 # Test for problem with broadcast dimensions
 @test lsim(sys, zeros(1, 5), 0:0.2:0.8)[1][:] == zeros(5)