h=0.01;
x(1,:)=[6400,0.7,8,0.32];
for k=1:1:101
t(k)=(k-1)*h;
xx=x(k,:)';
K1=ff(t(k),xx);
K2=ff(t(k)+h/2,xx+(h/2)*K1);
K3=ff(t(k)+h/2,xx+(h/2)*K2);
K4=ff(t(k)+h,xx+h*K3);
x(k+1,:)=x(k,:)+(h*(K1+2*K2+2*K3+K4)/6)';
end
for k=1:1:101
y(k,1)=x(k,1)*cos(x(k,2));
y(k,2)=x(k,1)*sin(x(k,2));
end
plot(y(:,1),y(:,2),'r-');
