2.2.m

来自「利用龙格库塔法（等步长」· M 代码 · 共 34 行

t0=0;
tend=1;
u0=[1,0,0]'
for j=2:5
    N=125*2^(j-1);
    h=(tend-t0)/N;
    result=[u0'];
    time=[t0];
    u=u0;
    t=t0;
    for k=1:N
        k1=fvdp1(t,u);
        k2=fvdp1(t+h,u+h*k1);
        k3=fvdp1(t+h/2,u+h*k1/4+h*k2/4);
        u=u+h*(k1+k2+4*k3)/6;
        t=t+h;
        result=[result;u'];
        time=[time;t];

    end
        plot(result(:,1))
        hold on
        plot(result(:,2))
        hold on
        plot(result(:,3))
        hold on
        loglog(time,result(:,1))
        hold on
        loglog(time,result(:,2))
        hold on
        loglog(time,result(:,3))
end