I have multiple ODE problems that I am solving. where I need the solutions (u) and derivative of the solutions (du). For smaller ODEs it is practical for me to do the following

```
using DifferentialEquations
function SB(du,u,p,t)
du[1]=@. u[2]
du[2]=@. ((-0.5*u[2]^2)*(3-u[2]/(p[4]))+(1+(1-3*p[7])*u[2]/p[4])*((p[6]-p[5])/p[2]+2*p[1]/(p[2]*p[9]))*(p[9]/u[1])^(3*p[7])-2*p[1]/(p[2]*u[1])-4*p[3]*u[2]/(p[2]*u[1])-(1+u[2]/p[4])*(p[6]-p[5]+p[10]*sin(2*pi*p[8]*t))/p[2]-p[10]*u[1]*cos(2*pi*p[8]*t)*2*pi*p[8]/(p[2]*p[4]))/((1-u[2]/p[4])*u[1]+4*p[3]/(p[2]*p[4]))
end
R0=2e-6
ps=250e3
f=2e6
u0=([R0 0])
tspan=(0,100/f)
p=[0.0725, 998, 1e-3,1481, 0, 1.01e5,7/5,f, R0, ps]
prob = ODEProblem(SB,u0,tspan,p)
@time u = solve(prob,Tsit5(),alg_hints=[:stiff],saveat=0.01/f,reltol=1e-8,abstol=1e-8)
t=u.t
u2=@. ((-0.5*u[2,:]^2)*(3-u[2,:]/(p[4]))+(1+(1-3*p[7])*u[2,:]/p[4])*((p[6]-p[5])/p[2]+2*p[1]/(p[2]*p[9]))*(p[9]/u[1,:])^(3*p[7])-2*p[1]/(p[2]*u[1,:])-4*p[3]*u[2,:]/(p[2]*u[1,:])-(1+u[2,:]/p[4])*(p[6]-p[5]+p[10]*sin(2*pi*p[8]*t))/p[2]-p[10]*u[1,:]*cos(2*pi*p[8]*t)*2*pi*p[8]/(p[2]*p[4]))/((1-u[2,:]/p[4])*u[1,:]+4*p[3]/(p[2]*p[4]))
```

where u2 is bascially du[2] in the SB function. This quickly becomes impractical as the size of my ODEs grow (>500 coupled ODEs with >500X500 matrices). Is there way to ask DifferentialEquations.jl package (or any other way) to export du[i]s as it is solving the ODEs? as a side note, it is important to mention that i need every time step that I am asking the solver to save the output at (saveat=0.01/f). I learned that using `integrator`

and `get_du!`

I should be able to do the job.

```
Rdot=zeros(50001,2)
u = init(prob,SSPRK22(),dt=1e-9,reltol=1e-8,abstol=1e-8)
solve!(u)
get_du!(Rdot,u)
```

However the obtained results from get_du! do not match u2. Am I missing something? appreciate the help.