matlab代码:
function z=fm(x,y)
z=x^2-y^2;
end
function [x,y]=RK(f,a,b,h,x0,y0)
n=ceil((b-a)/h);
x=zeros(n+1,1);
y=zeros(n+1,1);
x(1)=x0;
y(1)=y0;
for i=1:n
x(i+1)=x(i)+h;
k1=h*feval(f,x(i),y(i));
k2=h*feval(f,x(i)+h/2,y(i)+k1/2);
k3=h*feval(f,x(i)+h/2,y(i)+k2/2);
k4=h*feval(f,x(i)+h,y(i)+k3);
y(i+1)=y(i)+(k1+2*k2+2*k3+k4)/6;
end
end
运行结果: