2.1.m

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

global epsilon
epsilon=1;
t0=0;
tend=1;
u0=[1,0]'
for j=6:-1:1
    N=10*2^(j-1);
    h=(tend-t0)/N;
    result=[u0'];
    time=[t0];
    u=u0;
    t=t0;
    for k=1:N
        k1=fvdp(t,u);
        k2=fvdp(t+h,u+h*k1);
        k3=fvdp(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];
        if k==320
            result1=result(320,1);
        end
        
        if k==N
            result(:,1)=result(:,1)-result1;
            plot(time,result(:,1))
            hold on
            loglog(time,result(:,1))
        end
    end
end