% Integrates the system with RK4
% X0 : row vector initial condition = [X0,Y0,Z0]
% dt : scalar, time step for integration
% iter : scalar, number of iterations for integration
% XRK4 : matrix size (iter+1,3), contain [X(t),Y(t),Z(t)] for iter+1 values of t

function [XRK4] =RK4(X0,dt,iter)

% Initialising matrices
XRK4=zeros(iter+1,size(X0,2)); % Initialise matrices that will contain [X;Y;Z] vectors for all time steps
XRK4(1,:)=X0; % Setting initial state

for n=1:iter
	t=n*dt; 

    % RungeKutta 4
    k1= dt*diffEq(t,XRK4(n,:))';
    k2= dt*diffEq(t+dt/2.,XRK4(n,:)+ k1/2.)';
    k3= dt*diffEq(t+dt/2.,XRK4(n,:)+ k2/2.)';
    k4= dt*diffEq(t+dt,XRK4(n,:)+ k3)';
    XRK4(n+1,:)= XRK4(n,:)+(k1+2*k2+2*k3+k4)/6.;
end

end


