function [x,y] = RKs3o3B(fun,a,b,y0,N)
% Runge Kutta with 3 stages 3rd order
% num_ODEs_gr.pdf page 44 (second method)
% 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, y(n) + h*(-k1+2*k2)) ;
     y(n+1) = y(n) + h*(k1+4*k2+k3)/6 ;
end    

end

