function [x,Y]=RK2_sys3A(fun,Y0,x0,xN,N)
% Runge-Kutta method 2nd order 2 stages (s=2) for a system of 3 odes
% 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)x3 matrix that contains the values of Y(x) at the partition points
% comparing with RK2_sys3.m

a21=1; b1=1/2; b2 = 1/2;
% step length
h = (xN-x0)/N ;

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

% initial condition
Y(1,:) = Y0 ; % lines 27, 28, 29

for n=1:N
    x(n+1) = x(n) + h ;
    k1=feval(fun,Y(n,:)); % k1 is a vector
    k2=feval(fun,Y(n,:)+h*a21*k1); % k2 is a vector
    
    Y(n+1,:) = Y(n,:) + h*(b1*k1+b2*k2) ;
   end    
end       
