%Solution of First Order Coupled DEs
%-------------------------
%dx/dt=xy+t             equation I
%dy/dt=x-t              equation II
%at t=0, x=0 and y=1    initial conditions
%-------------------------

f1=@(t,x,y) x.*y+t;                   %  equation I
f2=@(t,x,y) x-t;    %  equation II
m=200;                              %number of x and t values
t=linspace(0.0,2,m);
x=zeros(1,m);
y=zeros(1,m);
x(1)=0;                             %initial conditions
y(1)=1;
h=t(2)-t(1);                        %interval

for n=1:m-1                         %4th order R-K Method
    k1x=h*f1(t(n),x(n),y(n));       %k values for equation I
    k1y=h*f2(t(n),x(n),y(n));       %k values for equation II
    k2x=h*f1(t(n)+h/2,x(n)+k1x/2,y(n)+k1y/2);
    k2y=h*f2(t(n)+h/2,x(n)+k1x/2,y(n)+k1y/2);
    k3x=h*f1(t(n)+h/2,x(n)+k2x/2,y(n)+k2y/2);
    k3y=h*f2(t(n)+h/2,x(n)+k2x/2,y(n)+k2y/2);
    k4x=h*f1(t(n)+h,x(n)+k3x,y(n)+k3y);
    k4y=h*f2(t(n)+h,x(n)+k3x,y(n)+k3y);

    x(n+1)=x(n)+(1/6)*(k1x+2*k2x+2*k3x+k4x);
    y(n+1)=y(n)+(1/6)*(k1y+2*k2y+2*k3y+k4y);
end;
subplot(1,2,1)
plot(t,x,t,y,'LineWidth',2);        % graphs t vs x and t vs y
legend('t vs x','t vs y','Location','North');
subplot(1,2,2);
plot(x,y,'LineWidth',2);            %graph x vs y
legend('x vs y','Location','North');
clear all;