Apologies if there has been a trivial mistake or otherwise in what follows.

MWE:

```
using Symbolics, Plots, DifferentialEquations, ModelingToolkit,
LaTeXStrings, Measures, DynamicalSystems
function forced_osc!(du, u, p, t) #the system
ω = p[1]; F = p[2];
du[1] = u[2]
du[2] = - sin(u[1])*ω^2 + F*cos(ω*t)*ω^2
return nothing
end
u0 = [0, 1e-4]
diffeq = (alg = Rodas4P(), abstol = 1e-8, reltol = 1e-8, maxiters = 1e8) #Solver specifications
forcedoscillator = ContinuousDynamicalSystem(forced_osc!, u0, [1,0.5]; diffeq)
ψ, t = trajectory(forcedoscillator6, 800; Δt = 0.1)
PS = scatter(ψ6[:, 1], ψ6[:, 2], lw = 0, markersize = 0.5, size = (700, 400), markerstrokewidth = 0, framestyle = :box, label = L"F = 0.5", legendposition = :topright)
```

Now, this runs perfectly and so and so. The issue is that the solution differs based on what algorithm I use.

Examples:

With `Rodas4P()`

:

With `RK4()`

:

With `Tsit5()`

(default):

With `AutoVern7(Rodas4())`

(recommended in the docs for lower tolerances):

All of these have differing long-time behaviours - is this expected instability associated with the solvers? I started with a stiff solver because the stiffness was unknown (and I had used this `Rodas4P()`

combination in a previous problem) - as far as I know, stiff solvers can be used for non-stiff problems but the efficiency drops. Since my system wasn’t taking much time, I didn’t reflect on it much.

Is there a reason/solution for this? I noted that the `Rodas4P()`

solution goes more towards the rest if the tolerances are further lowered to `1e-9`

. Since the rest of the solvers branch towards the right, I am guessing the `Rodas4P()`

solution is the least credible here?

I have also verified that this same thing also happens if one only uses these methods in `DifferentialEquations.jl`

only (i.e., without using `DynamicalSystems.jl`

).

Any help will be appreciated!

**Note:** It does seem like an inbuilt solver error - lowering the tolerances to `1e-9`

gives proper match between these solutions for around `tspan = (0, 850)`

.