Computing the adjoint for an equation integrated backward in time

[guachoperez]

I want to integrate a system of ODEs backward in time over some interval. Then I want to measure the initial boundary conditions and differentiate with respect to the parameters of the system. I think that the following code snippet is a very basic illustration of what I want to achieve:

f(u, p, t) = p[1] * u
sol(p) = solve(ODEProblem(f, [0.1], (1.0, 0.0), p), dt = -1)
gradient(p -> sum(sol(p)), [1.01])

The final gradient is just an example of what I want to do, but both this demo and my actual application run into the same problem. Namely, I get the warning:

Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.

Is there a way to achieve what I want to do without having to manually subsitute u = 1 - t into the system?