multistep.m

来自「various matlab programs to slove various」· M 代码 · 共 23 行

function [t,u] = multistep (a,b,tf,t0,u0,h,fun,dfun,tol,itmax)
y = u0;  t = t0;  f = eval (fun); p = length(a) - 1; u = u0;
nt = fix((tf - t0 (1) )/h);
for k = 1:nt
    lu=length(u);
    G = a' *u (lu:-1:lu-p)+ h * b(2:p+2)' * f(lu:-1:lu-p);
    lt = length(t0); t0 = [t0; t0(lt)+h]; unew = u (lu); 
    t = t0 (lt+1); err = tol + 1; it = 0;
    while (err > tol) & (it <= itmax)
      y = unew; den = 1 - h * b (1) * eval(dfun);
      fnew = eval (fun);
      if den == 0
         it = itmax + 1;
      else
         it = it + 1;
         unew = unew - (unew - G - h * b (1) * fnew)/den;
         err = abs (unew - y);
      end
    end
    u = [u; unew]; f = [f; fnew];
end   
t = t0;