2.5.m

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

global ka;
t0=0;
tend=1.5;
A=[0.020,0.065];
alpha=[pi/6,pi/4,pi/3];
for m=1:2
    ka=A(m);
    for n=1:3        
        u0=[0,1.5,20*cos(alpha(n)),20*sin(alpha(n))]';
        N=100;
        h=(tend-t0)/N;
        result=[u0'];
        time=[t0];
        u=u0;
        t=t0;   
        for k=1:N
            k1=fvdp3(t,u);
            k2=fvdp3(t+h,u+h*k1);
            k3=fvdp3(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),result(:,2))
            hold on
    end
end