Today I found that the solver Rodas4P2 is working best for my model. For an tolerance of 1e-7 it is 4 times faster than Simulink, operating at a tolerance of 1e-3.

But what is the difference between Rodas4, Rodas4P and Rodas4P2?

I know that Rodas is described in literature, but I could not find any paper explaining Rodas4P or Rodas4P2… Is there any paper about it?

If you’re using a step input, you should probably be using tstops or a ContinuousCallback to tell the solver where the discontinuity is (see Handling Instability When Solving ODE Problems).

Just tested to add tstops. It works, but it is not sufficient to add a tstop at the time of the discontinuity, but I need one more for a perfect result. Example (assuming the step happens at 200s):

function simulate(sys, ω, ω0, U, U0, D; p = [ U => U0 ], prn=true, duration=1000)
u0 = [ω => ω0, # this is a fixed initial value
D(ω) => 0.0] # the value for the derivative is an inital estimate
# simulation parameters
tspan = (0.0, duration)
dt = 0.01
ts = 0:dt:duration
# run the simulation
prob = ODEProblem(sys, u0, tspan, p)
sol = solve(prob, Rodas5(), tstops=[200, 200+dt], saveat=ts)
end