龙格-库塔

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

运行结果：