fvdp3.m

利用龙格库塔法（等步长

function rhs=fvdp3(t,u)
global ka
rhs=[u(3)
    u(4)
    -ka*u(3)*sqrt(u(3)^2+u(4)^2)
    -9.81-ka*abs(u(4))*sqrt(u(3)^2+u(4)^2)];