function [x,y] = RKs4(fun,A,b,c,xstart,xend,y0,N)
% Runge Kutta with 4 stages 4th order 
% A, b ,c Butcher array
% num_ODEs_gr.pdf page 45
% fun the right hand side function f (ODE y'=f(x,y))
% Integration interval [xstart,xend]
% Initial condition y0
% Number of (partition) points N+1

% x a column vector contains the partition points
% y a column vector contains the numerical solution at the partition points
% initialize with zeros
x = zeros(N+1,1) ;
y = zeros(N+1,1) ;

% N subintervals
% find the step length
h = (xend-xstart)/N ;

x(1) = xstart ; % integration interval 
y(1) = y0 ; % initial condition [xstart,xend]
for n = 1:N
     x(n+1) = x(n) + h ;
     k1 = feval(fun, x(n), y(n)) ;
     k2 = feval(fun, x(n)+ c(2)*h, y(n) + h*A(2,1)*k1) ;
     k3 = feval(fun, x(n)+ c(3)*h, y(n) + h*(A(3,1)*k1+A(3,2)*k2) ;
     k4 = feval(fun, x(n)+ c(4)*h, y(n) + h*(A(4,1)*k1+ A(4,2)*k2+A(4,3)*k3)) ;
     y(n+1) = y(n) + h*(b1*k1+b2*k2+b3*k3+b4*k4) ;
end    
end

