f1=@(t,x,y) x.*y+t;
f2=@(t,x,y) x-t;
m=200;
t=linspace(0.0,2,m);
x=zeros(1,m);
y=zeros(1,m);
x(1)=0;
y(1)=1;
h=t(2)-t(1);
for n=1:m-1
k1x=h*f1(t(n),x(n),y(n));
k1y=h*f2(t(n),x(n),y(n));
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);
legend('t vs x','t vs y','Location','North');
subplot(1,2,2);
plot(x,y,'LineWidth',2);
legend('x vs y','Location','North');
clear all;