function [x,y1,y2,y3]=RK2_sys3(fun,y10,y20,y30,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
% y10,y20,y30 the initial conditions y1(x0), y2(x0), y3(x0)
% [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 
% y1 a vector that contains the values of y_1(x) at the partition points
% y2 a vector that contains the values of y_2(x) at the partition points
% y3 a vector that contains the values of y_3(x) at the partition points

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);
y1 = zeros(N+1,1);
y2 = zeros(N+1,1);
y3 = zeros(N+1,1);
% starting point of the integration interval
x(1) = x0 ;

% initial condition
y1(1) = y10 ;
y2(1) = y20;
y3(1) = y30;

for n=1:N
    x(n+1) = x(n) + h ;
    [k1,l1,m1] = feval(fun,y1(n),y2(n),y3(n));
    [k2,l2,m2] = feval(fun,y1(n)+h*a21*k1,y2(n)+h*a21*l1,y3(n)+h*a21*m1);
    
    y1(n+1) = y1(n) + h*(b1*k1+b2*k2) ;
    y2(n+1) = y2(n) + h*(b1*l1+b2*l2) ;
    y3(n+1) = y3(n) + h*(b1*m1+b2*m2) ;
 end    
end       
