function ode_msd()

clf

%%%%%%%%%%%%%%%%%%%%
%%% Euler's method
%%%%%%%%%%%%%%%%%%%

%% initial conditions
x0 = 1;
v0 = 0;
t = 0;

%% pack the data into a vector
y = [x0;
     v0];

%% integration parameters
dt = 0.05;


for i = 1:1000
  dy_dt = deriv(t,y);
  y = y + dy_dt*dt;  
  t = t+dt;
  plot(t,y(1),'.');
  hold on
end



%%%%%%%%%%%%%%%%%%%%%%%%
%%% Midpoint method
%%%%%%%%%%%%%%%%%%%%%%%%
%% initial conditions
x0 = 1;
v0 = 0;
t = 0;

%% pack the data into a vector
y = [x0;
     v0];


for i = 1:1000
  dy_dt = deriv(t,y);
  yH    = y + dy_dt*dt/2;
  dy_dt = deriv(t,yH);
  y     = y + dy_dt*dt;
  
  t = t+dt;
  plot(t,y(1),'r.');
  hold on
end






%%%%%%%%%%%%%%%%%%%%%%%%
%%% ODE45
%%%%%%%%%%%%%%%%%%%%%%%%
%% initial conditions
x0 = 1;
v0 = 0;

%% pack the data into a vector
y = [x0;
     v0];

options = odeset('AbsTol',1e-6,'RelTol',1e-6);
[T,Y] = ode45(@deriv,[0 50],y,options);


plot(T,Y(:,1),'k','LineWidth',2)


set(gca,'FontSize',14);
xlabel('time')
ylabel('position')
grid on







%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dy_dt = deriv(t,y)

  %% unpack vector
  x = y(1);
  v = y(2);
  
  dx_dt = v;
  dv_dt = -x;  %% spring equation F=ma=-kx
  
  dy_dt = [dx_dt;
	   dv_dt];

