Here is a (mildly) interesting example x''+x'+x=\pm p_1 where the sign of p_1 changes when a switching manifold is encountered at x=p_2. To make things more interesting, consider hysteresis in the switching manifold such that p_2\mapsto -p_2 whenever the switching manifold is crossed.

The code is relatively straightforward; the StaticArrays/SVector/MVector can be ignored, they are only for speed.

```
using OrdinaryDiffEq
using StaticArrays
f(x, p, t) = SVector(x[2], -x[2]-x[1]+p[1]) # x'' + x' + x = ±p₁
h(u, t, integrator) = u[1]-integrator.p[2] # switching surface x = ±p₂;
g(integrator) = (integrator.p .= -integrator.p) # impact map (p₁, p₂) = -(p₁, p₂)
prob = ODEProblem(f, # RHS
SVector(0.0, 1.0), # initial value
(0.0, 100.0), # time interval
MVector(1.0, 1.0)) # parameters
cb = ContinuousCallback(h, g)
sol = solve(prob, Vern6(), callback=cb, dtmax=0.1)
```

Then plot `sol[2,:]`

against `sol[1,:]`

to see the phase plane - a nice non-smooth limit cycle in this case.

Note that if you try to use interpolation of the resulting solution (i.e., `sol(t)`

) you need to be very careful around the points that have a discontinuous derivative as the interpolant goes a little awry. That’s why I’ve used `dtmax=0.1`

to get a smoother solution output in this case. (I’m probably not using the most appropriate integrator either but it’s the one that I was using in a previous piece of code that I copied-and-pasted )

If you are interested in the dynamics of piecewise differential equations there are plenty of gotchas to watch out for, particularly if your system is a Filippov system; effects like sliding require a bit more care than shown in the code above. (E.g., I’m not sure if DifferentialEquations.jl can handle a system that changes from an ODE to a DAE part way through solving, such as happens with the sliding vector fields that come from Filippov systems; ever tried it @ChrisRackauckas?)