2.1.1.m

利用龙格库塔法（等步长

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 j==6
            result1=result(320);
        end
        result=result-result(320);
        plot(loglog(t,result(:,1)),f9,fdsdf)
    end