function [t,y] = rk2(f,t0,tend,y0,N)
%
% [tn,yn] = rk2(f,t0,tend,y0,N)
%
% Solves the ODE
%
% y'=f(t,y), y(t0)=y0, t0<= t <= tend
%
% by using N constant stepsizes, using 
% a second order RK method. 
%
m = length(y0);

h = (tend - t0)/N;  % Stepsize

t = zeros(1,N+1);       % Store the timesteps 
y = zeros(m,N+1);       % Store the numerical solutions

t(1) = t0; 
y(:,1) = y0;

% The main loop (time stepping)
for n=1:N
    t(n+1) = t(n)+h; 
    k1 = f(t(n),y(:,n)); 
    k2 = f(t(n)+h,y(:,n)+h*k1);
    y(:,n+1)=y(:,n)+0.5*h*(k1+k2);
end




