f1=@(t,x1,x2) x2;
f2=@(t,x1,x2) -0.1*x2.^2-0.6*x1;
m=200;
t=linspace(0.0,20,m);
x1=zeros(1,m);
x2=zeros(1,m);
x1(1)=1;
x2(1)=0;
h=t(2)-t(1);
for n=1:m-1
k11=h*f1(t(n),x1(n),x2(n));
k12=h*f1(t(n)+h/2,x1(n)+k11/2,x2(n));
k13=h*f1(t(n)+h/2,x1(n)+k12/2,x2(n));
k14=h*f1(t(n)+h,x1(n)+k13,x2(n));
x1(n+1)=x1(n)+(1/6)*(k11+2*k12+2*k13+k14);
k21=h*f2(t(n),x1(n),x2(n));
k22=h*f2(t(n)+h/2,x1(n),x2(n)+k21/2);
k23=h*f2(t(n)+h/2,x1(n),x2(n)+k22/2);
k24=h*f2(t(n)+h,x1(n),x2(n)+k23);
x2(n+1)=x2(n)+(1/6)*(k21+2*k22+2*k23+k24);
end;
subplot(1,2,1)
plot(t,x1,t,x2,'LineWidth',2);
legend('t vs x','t vs dx/dt','Location','SouthWest');
subplot(1,2,2);
plot(x1,x2,'LineWidth',2);
legend('X1 vs X2','Location','SouthWest');
clear all;