nark4v.m

自己编写的变步长龙格库塔方法解强非线性微分方程

function [t,y]=nark4v(tspan,y0,e,h)
if nargin<4
    h=(tspan(2)-tspan(1))/100;
end
n=1;
t(n)=tspan(1);
y(:,n)=y0;
[y1,y2]=comput(t(n),y(:,n),h);
while t(n)<tspan(2)-eps
    if abs(y2(2)-y1(2))/10>e
        while (abs(y2(2)-y1(2))+abs(y2(1)-y1(1)))/10>e
            h=h/2;
            [y1,y2]=comput(t(n),y(:,n),h);
        end
    else
        while (abs(y2(2)-y1(2))+abs(y2(1)-y1(1)))/10<=e
            h=2*h;
            [y1,y2]=comput(t(n),y(:,n),h);
        end
        h=h/2;
        h=min(h,tspan(2)-t(n));
        [y1,y2]=comput(t(n),y(:,n),h);
    end
    n=n+1;
    t(n)=t(n-1)+h;
    y(:,n)=y2;
    [y1,y2]=comput(t(n),y(:,n),h);
end
t=t'; y=y';