function [x,Y]=RK4_sys(fun,Y0,x0,xN,N,A,b,c)
% Runge-Kutta method 4 stages (s=4) for any system of  odes
% A (sxs), b (1xs), c (1xs) the coefficients of the method must be given by the user
% Input arguments 
% fun the right hand side function
% Y0 the initial condition vector
% [x0,xN] integration interval
% N the number of subintervals (partition of N+1 points)
% Output arguments
% x a vector with the points of the partition x_0,x_1,x_2,...,x_N 
% Y is (N+1)xM matrix that contains the values of Y(x) at the partition points
% M is the length of the initial condition vector (Y0)

M = length(Y0) ;

% step length
h = (xN-x0)/N ;

% creates the vectors to be used
x = zeros(N+1,1);
Y = zeros(N+1,M); 
% starting point of the integration interval
x(1) = x0 ;

% initial condition
Y(1,:) = Y0 ; 

for n=1:N
    x(n+1) = x(n) + h ;
    k1=feval(fun,Y(n,:));
    k2=feval(fun,Y(n,:)+h*A(2,1)*k1);
    k3=feval(fun,Y(n,:)+h*(A(3,1)*k1+A(3,2)*k2)); 
    k4=feval(fun,Y(n,:)+h*(A(4,1)*k1+A(4,2)*k2+A(4,3)*k3));
    Y(n+1,:) = Y(n,:) + h*(b(1)*k1+b(2)*k2+b(3)*k3+b(4)*k4) ;
   end    
end       
