function [tb,x,y] = RK4System(int,x0,y0,nsteps)
% Runge-Kutta System
%Input: Interval [ta,T], initial condition x0,  and number of steps n
%Output: the approximate solution x, at time T
% Example usage: [T,x,y]=RK4System([0 1],1,0,100);


ta = int(1); tb = int(2);
h = (tb - ta)/nsteps; %%step size h=(T-ta)/n
time=zeros(nsteps,1);
x(1) = x0;
y(1) = y0;
time(1) = ta;


format long 

for j = 1:nsteps
    K11 = h*xdot(time(j),x(j),y(j));
    K12 = h*ydot(time(j),x(j),y(j));
    
    K21 = h*xdot(time(j)+h/2,x(j)+K11/2,y(j)+K12/2);
    K22 = h*ydot(time(j)+h/2,x(j)+K11/2,y(j)+K12/2);
    
    K31 = h*xdot(time(j)+h/2,x(j)+K21/2,y(j)+K22/2);
    K32 = h*ydot(time(j)+h/2,x(j)+K21/2,y(j)+K22/2);
    
    K41 = h*xdot(time(j)+h,x(j)+K31,y(j)+K32);
    K42 = h*ydot(time(j)+h,x(j)+K31,y(j)+K32);
    
    x(j+1) = x(j) + (K11 + 2*K21 + 2*K31 + K41)./6;
    y(j+1) = y(j) + (K12 + 2*K22 + 2*K32 + K42)./6;
    
    time(j+1) = time(j) + h;
end

figure(1)
    title('System of 2 equations with RK4')
    hold all
    plot(time,x,'.-.','LineWidth',1,...
                'MarkerEdgeColor','k',...
                'MarkerFaceColor','g',...
                'MarkerSize',10);

    plot(time,y,'-','LineWidth',2,...
                'MarkerEdgeColor','k',...
                'MarkerFaceColor','g',...
                'MarkerSize',10);
legend('x','y')

figure(2)
hold all
plot(x,y,'o')
axis on
title('Phase plane x-y')


end



function z1 = xdot(t,x,y) %the RHS of the ODE
z1 = x - y + t*(2-t*(1+t));
end

function z2 = ydot(t,x,y)
z2 = x + y - t^2*(4-t);
end
