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
x = x0;
y = y0;
t = ta;
format long 

for j = 1:nsteps
    K11 = h*xdot(t,x,y);
    K12 = h*ydot(t,x,y);
    
    K21 = h*xdot(t+h/2,x+K11/2,y+K12/2);
    K22 = h*ydot(t+h/2,x+K11/2,y+K12/2);
    
    K31 = h*xdot(t+h/2,x+K21/2,y+K22/2);
    K32 = h*ydot(t+h/2,x+K21/2,y+K22/2);
    
    K41 = h*xdot(t+h,x+K31,y+K32);
    K42 = h*ydot(t+h,x+K31,y+K32);
    
    x = x + (K11 + 2*K21 + 2*K31 + K41)./6;
    y = y + (K12 + 2*K22 + 2*K32 + K42)./6;
    
    t = t + h;
figure(1)
title('System of 2 equations with RK4')
plot(t,x,'.','LineWidth',1,...
                'MarkerEdgeColor','k',...
                'MarkerFaceColor','g',...
                'MarkerSize',10);

plot(t,y,'o','LineWidth',1,...
                'MarkerEdgeColor','k',...
                'MarkerFaceColor','g',...
                'MarkerSize',10);

hold all
end



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

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

