ODE = @(t,x)[x(2);x(1)];
t = linspace(0,1,1000);
% t = [0 5];
[t,x] = ode23(@(t,x)[x(2);x(1)],t,[1;0]);
[ts,xs] = ode23(@(t,x)[x(2);x(1)],[0,1],[1;0]);
%[t,x] = ode45(@(t,x)[x(2);x(1)],t,[1;0]);
%[ts,xs] = ode45(@(t,x)[x(2);x(1)],[0,1],[1;0]);
vs_interp = interp1(ts,xs(:,2),t);
xs_interp = interp1(ts,xs(:,1),t);


figure(1)
plot(t,x)
xlabel('time'); ylabel('position and speed')
legend('position','speed')
legend('position','speed','location','northwest')
figure(2)
t_disp = t(1:end-1)+diff(t)/2;
a = diff(x(:,2))./diff(t);
ts_disp = ts(1:end-1)+diff(ts)/2;
as = diff(xs(:,2))./diff(ts);
plot(t_disp,a,ts_disp,as,'+-')
legend('small time steps','default time steps')
xlabel('time'); ylabel('acceleration')

if 0
  figure(3)
  plot(t,x(:,2),ts,xs(:,2),'+-')
  xlabel('time'); ylabel('speed')

  figure(4)
  plot(t,x(:,1)-xs_interp)
endif