impeuler.m

来自「matlab codes for solving ODEs, including」· M 代码 · 共 16 行

function [tout,yout]=impeuler(f,tspan,y0,N)

h=(tspan(2)-tspan(1))/N;
t = tspan(1);
tout = t;
y=y0(:);
yout = y.';
for n=1:N
    f1 = feval(f,t,y);
    f2 = feval(f,t+h,y+h*f1);
    y = y+h*(f1+f2)/2;
    t = t+h;
    yout=[yout;y.'];
    tout=[tout;t];
end