```julia using OrdinaryDiffEqTsit5, ModelingToolkit, StaticArrays using ModelingToolkit: t_nounits as t, D_nounits as D @parameters σ ρ β @variables x(t) y(t) z(t) eqs = [D(D(x)) ~ σ * (y - x), D(y) ~ x * (ρ - z) - y, D(z) ~ x * y - β * z] @named sys = ODESystem(eqs, t) sys = structural_simplify(sys; split=false) u0 = SA[D(x) => 2f0, x => 1f0, y => 0f0, z => 0f0] p = SA[σ => 28f0, ρ => 10f0, β => 8f0 / 3f0] tspan = (0f0, 100f0) prob = ODEProblem{false}(sys, u0, tspan, p) sol = solve(prob, Tsit5()) ```