f=@(x,t)x^2*cos(t);
h=1/16;
x(1)=1/2;
t(1)=0;
T=8;
for n=1:T/h
    f1=f(x(n),t(n));
    f2=f(x(n)+0.5*h*f1,t(n)+0.5*h);
    f3=f(x(n)+0.5*h*f2,t(n)+0.5*h);
    f4=f(x(n)+h*f3,t(n)+h);
    x(n+1)=x(n)+(h/6)*(f1+2*f2+2*f3+f4);
    t(n+1)=t(1)+h*n;
end;
plot(t,x);