function [x,y] = RKs4o4C(fun,a,b,y0,N)
% Runge Kutta with 4 stages 4th order 
% num_ODEs_gr.pdf page 46 (3rd method)
% classical Runge-Kutta
% Integration interval [a,b]
% Initial condition y0
% Number of (partition) points N+1

x = zeros(N+1,1) ;
y = zeros(N+1,1) ;

h = (b-a)/N ;
x(1) = a ;
y(1) = y0 ;
for n = 1:N
     x(n+1) = x(n) + h ;
     k1 = feval(fun, x(n), y(n)) ;
     k2 = feval(fun, x(n)+h/2, y(n) + h*k1/2) ;
     k3 = feval(fun, x(n)+h/2, y(n) + h*k2/2) ;
     k4 = feval(fun, x(n)+h, y(n) + h*k3) ;
     y(n+1) = y(n) + h*(k1+2*k2+2*k3+k4)/6 ;
end    
end

