euler_r.m

matlab具有强大的数值计算能力

function [a,b]=Euler_r(ydot_fun,a0,b0,h,N)
%用改进Euler公式，其中
%ydot_fun为一阶微分方程的函数；
%x0,y0为初始条件
%h为区间长度
%N为区间的个数
%x为Xn构成的向量
%y为Yn构成的向量
a=zeros(1,N+1);b=zeros(1,N+1);
a(1)=a0;y(1)=b0;
for n=1:N
    a(n+1)=a(n)+h;
    bbar=b(n)+h*feval(ydot_fun,a(n),b(n));
    b(n+1)=b(n)+h/2 * (feval(ydot_fun,a(n),b(n))+feval(ydot_fun,a(n+1),bbar));
end