%4階經(jīng)典RK
function [x,y]=rgkl();
input('請(qǐng)按下列步驟輸入初始值!');
x0=input('x0:');
y0=input('y0:');
h=input('步長(zhǎng):');
e=input('最大范圍:=');

i=1:h:e;
n=max(size(i));
x(1)=x0;
y(1)=y0;
for t=1:n-1
k1=x(t)-4*y(t)*sqrt(abs(x(t)*y(t)));
k2=(x(t)+h/2)-4*(y(t)+h/2*k1)*sqrt(abs((x(t)+h/2)*(y(t)+h/2*k1)));
k3=(x(t)+h/2)-4*(y(t)+h/2*k2)*sqrt(abs((x(t)+h/2)*(y(t)+h/2*k2)));
k4=(x(t)+h)-4*(y(t)+h*k3)*sqrt(abs((x(t)+h)*(y(t)+h*k3)));
l1=-y(t)+4*x(t)*sqrt(abs(x(t)*y(t)));
l2=(-y(t)+h/2*l1)+4*(x(t)+h/2)*sqrt(abs((x(t)+h/2)*(y(t)+h/2*l1)));
l3=(-y(t)+h/2*l2)+4*(x(t)+h/2)*sqrt(abs((x(t)+h/2)*(y(t)+h/2*l2)));
l4=(-y(t)+h*l3)+4*(x(t)+h)*sqrt(abs((x(t)+h)*(y(t)+h*l3)));
x(t+1)=x(t)+h/6*(k1+2*k2+2*k3+k4);
y(t+1)=y(t)+h/6*(l1+2*l2+2*l3+l4);
end
plot(x,y)