euler1.m

有关matlab r2007的实用教程

function [xout,yout]=Euler1(fun,x0,xn,y0,n,eps)
%用欧拉方法计算常微分方程的初值问题
%x0,y0为初值条件
%xn,yn为x取值区间的最后一个节点的横坐标和纵坐标
% n为区间的等分数目
pow=1/3;
if nargin<5
    eps=1e-3;
end
x=x0;
hmax=(xn-x0)/n;
h=hmax/8;
y=y0(:);
chunk=128;
xout=zeros(chunk,1);
yout=zeros(chunk,length(y));
k=1;
xout(k)=1;
yout(k,:)=y.';
while (x<xn)&(x+h>x)
    if x+h>xn
        h=xn-x;
    end
    f=feval(fun,x,y);
    f=f(:);
    delta=norm(h*f,'inf');
    tau=eps*max(norm(y,'inf'),1.0);
    if delta<=tau
        x=x+h;
        y=y+h*f;
        k=k+1;
        if k>length(xout)
            xout=[xout; zeros(chunk,1)];
            yout=[yout; zeros(chunk,length(y))];
        end
        xout(k)=x;
        yout(k,:)=y.';
    end
    if delta~=0
        h=min(hmax,0.9*h*(tau/delta)^pow);
    end
end
xout=xout(k);
yout=yout(k);