How to use a callback in DifferentialEquations.jl to simulate a hammer blow to a mass in a spring-mass system?

I'm teaching Differential Equations and am using Pluto notebooks with DifferentialEquations.jl to supplement the course. We are currently trying to model a spring-mass system that is subject to an impulse of 10 Newton-seconds (hammer blow) at a specific instant of time. I believe I need to use a callback to simulate this, but don't know how to implement it.

The differential equation for this system is: u''(t) + 2u'(t) + 10u(t) = 10*delta(t-1), where 10*delta(t-1) is the delta function for the 10 Newton-second impulse at 1 second. The initial conditions are u(0) = u'(0) = 0.

What I've done so far is:

function hammer_time(du, u, p, t)
    -2*du - 10*u
end

begin
    du0 = 0.0
    u0 = 0.0
    tspan = (0.0, 10.0)
end

prob = SecondOrderODEProblem(hammer_time, du0, u0, tspan)

begin
    hit_times = [1.0]
    affect!(integrator) = integrator.u += 0.5

    cb = PresetTimeCallback(hit_times, affect!)
    sol = solve(prob, Tsit5(), callback = cb)
end

In the last code block I know the second line is wrong, but I don't know what it should be to model the hammer blow impulse. I assume it changes the energy state of the system but I don't know what to do to model that.

As far as I can tell, SecondOrderODEProblem does not play nicely with callbacks mutating states - it originates from some details under the hood when constructing this nice interface, not important for your use case though.

However, we can certainly avoid this altogether if we write the 2nd order ODE as a system of first order ODEs:

Let u = [x; v], so that u' = [v; -2v-10x]

Then we can define our system using the more common ODEProblem interface

using DifferentialEquations, Plots;

#x = u[1], v = u[2]
function hammer_time(u, p, t)
    return [u[2];
            -2*u[2] - 10*u[1]];
end

u0 = [0.0;
      0.0];
tspan = (0.0, 10.0);

condition(u, t, integrator) = (t == 1);
affect!(integrator) = (integrator.u[2] += 10);
cb = DiscreteCallback(condition, affect!);
prob = ODEProblem(hammer_time, u0, tspan)
sol = solve(prob, Tsit5(), callback=cb, tstops=[1.0])
plot(sol, labels=["x" "v"])

This has the effect that at time 1, we increase the velocity (u[2]) by 10 units.

The above code yields plot of position and velocity as a function of time.

I'm sure you already found the documentation for callbacks, but here it is again for posterity. https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/#Using-Callbacks