naeuler2s.m

Matlab数值解算法实现代码

function [x,y]=naeuler2s(dyfun,xspan,y0,h)
%用途: 2阶改进euler法解常微分方程组y'=f(x,y),y(x0)=y0
%格式: [x,y]=naeuler2s(dyfun,xspan,y0,h),dyfun为函数f(x,y)，xspan为求解区间[x0,xN]
%      y0为初值向量y(x0)，h为步长,x返回节点，y返回数值解
x=xspan(1):h:xspan(2);
y=zeros(length(y0),length(x));
y(:,1)=y0(:);
for n=1:length(x)-1
    k1=feval(dyfun,x(n),y(:,n));
    y(:,n+1)=y(:,n)+h*k1;
    k2=feval(dyfun,x(n+1),y(:,n+1));
    y(:,n+1)=y(:,n)+h*(k1+k2)/2;
end
x=x';y=y';