DifferentialEquations.jl has an expressive callback system which allows for customizable transformations of the solver behavior. DiffEqCallbacks.jl is a library of pre-built callbacks which makes it easy to transform the solver into a domain-specific simulation tool.

For information on using the package, see the stable documentation. Use the in-development documentation for the version of the documentation, which contains the unreleased features.

Here we solve the harmonic oscillator:

u0 = ones(2)
function f(du, u, p, t)
    du[1] = u[2]
    du[2] = -u[1]
end
prob = ODEProblem(f, u0, (0.0, 100.0))

However, this problem is supposed to conserve energy, and thus we define our manifold to conserve the sum of squares:

function g(resid, u, p, t)
    resid[1] = u[2]^2 + u[1]^2 - 2
    resid[2] = 0
end

To build the callback, we just call

cb = ManifoldProjection(g)

Using this callback, the Runge-Kutta method Vern7 conserves energy. Note that the standard saving occurs after the step and before the callback, and thus we set save_everystep=false to turn off all standard saving and let the callback save after the projection is applied.

sol = solve(prob, Vern7(), save_everystep = false, callback = cb)
@test sol[end][1]^2 + sol[end][2]^2 ≈ 2