No, second order is not really all that high. You may still experience numerical difficulties due to poor scaling though, even for this simple system. Consider what happens if you increase `w1`

, for `w1=1e6`

, the solver exits with

```
┌ Warning: dt(8.881784197001252e-16) <= dtmin(8.881784197001252e-16) at t=1.9999999999999998, and step error estimate = 8.930759781475992. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase ~/.julia/packages/SciMLBase/s9wrq/src/integrator_interface.jl:599
retcode: DtLessThanMin
```

Already for `w1=1000`

, you can see part of the problem in the solution:

```
model, u, y = model_bessel(2, 100.0, w1=1e3)
prob = ODEProblem(model, [], tspan)
sol = solve(prob, Tsit5())
plot(sol)
```

The derivative state becomes enormous, much larger that the position state.

Compare this to a balanced representation:

```
using ControlSystemsMTK, ControlSystemsBase
function balanced_model_bessel(step_time, step_size; w1=3.14, y0=0.0, yd0=0.0)
println("Defining the model...")
@variables t u(t)=0 y(t)=y0 yd(t)=yd0
@parameters w=w1
H = tf(w1^2, [0.618, 1.3617w1, w1^2])
Hss = ss(H) # This performs balancing
@named filter = ODESystem(Hss)
eqs = [filter.input.u ~ u
u ~ ifelse(t > step_time, step_size, 0.0)]
@named sys = ODESystem(eqs, t, systems=[filter])
println("Model defined...")
sys = structural_simplify(sys)
println("Model simplified!")
sys, u
end
model, u = balanced_model_bessel(2, 100.0, w1=1e3)
prob = ODEProblem(model, [], tspan)
sol = solve(prob, Tsit5())
@show length(sol.t)
plot(sol)
```

Here, the two state variables are of comparable magnitude, making the job for the solver much easier.

To plot the correctly scaled filter output in this case, you’d have to plot `plot(sol, idxs=model.filter₊output₊u)`

, since the output is scaled w.r.t. the state.