Second derivative of an ODE

I have an initial value problem that I solve with Matlab's ode solver: dy/dx=f(x,y) with f(x_0)=y_0, where y is in R^2 and x is the time dimension. Apart from the solution I am also interested in the second derivative d^2y/dx^2. From what I understand of ODE-solvers, most solvers need this second derivative to find a good approximation of the initial value problem. However, as far as I can see the ODE-solvers do not provide me with the possibility to generate these series as an output. Obviously I could use the solution of the initial value problem to approximate the second derivative, but this will likely be less accurate than Matlab's approximation within the ODE-solver. Does anybody know how I can get the second derivative as an output?

Answer 1

[Grzegorz Knor]

Perhaps this simple source code will help you solve your problem:

function dx = ode_test(t,x)
dx(1) = sin(t);
dx(2) = x(1); 
dx(3) = x(2); 
dx = dx(:) ; 

[T,Y] = ode45(@ode_test,[0 2*pi],[-1 0 1])
plot(T,Y)