我需要修复代码以利用$2$阶段多步法：

$$y_{n+2} - y_n = h\left[(1/3)f_{n+2} + (4/3)f_{n+1} + (1/3)f_n\right]$$

由于这是一种隐式方法，因此我使用 Newton-Raphson 方法来最终确定$y_{n+2}$.

下面，是我尝试在第一个近似步骤中实现已经给定的 rk4 方法（runge-kutta 4）的代码：

function [tout, yout] = newcorrect(FunFcn,t0,tfinal,step,y0)
maxiter = 1000;
tolnr   = 1e-9;
diffdelta = 1e-6;
stages=2;
[tout,yout]=rk4(FunFcn,t0,t0+(stages-1)*step,step,y0);
tout=tout(1:stages);
yout=yout(1:stages);
t = tout(stages);
y = yout(stages).';
% The main loop
 while abs(t- tfinal)> 1e-6 
    if t + step > tfinal, step = tfinal - t; end
    t = t + step;
    yn0 = y;
    ynf = yn0;
    yn = inf;
    iter = 0;
    while (abs(yn - ynf)>= tolnr) && (iter < maxiter)
        df = 1/diffdelta * (feval(FunFcn,t, yn0+diffdelta) - feval(FunFcn, t, yn0));
        yn = yn0 - 1/(1/3*step*df - 1) * (4/3*step*feval(FunFcn,tout(end),yout(end)) + 1/3*step*feval(FunFcn,tout(end-1),yout(end-1)) + 1/3*step*feval(FunFcn, t, yn0) -yn0 + yout(end-1));
        ynf = yn0;
        yn0 = yn;
        iter = iter + 1;
    end
    y = yn;
    tout = [tout; t];
    yout = [yout; y.'];
 end
end

这用于通过实验证明，当您将步长除以 2 时，每个连续不同步长的最大绝对误差的分数$h$大约等于$2^{-p}$在哪里$p$是方法的顺序。该方法被证明是有序的$4$. 该练习希望最初的起点是$k0=2$因此我创建了以下脚本：

clear all;
t0 = 1;
tfinal = 3;
y1 = 2;
tout  =t0:0.01:tfinal;
k0 = 2;
kf = input('enter final k:')
for k = k0:kf
  h(k-1) = 2^(-k);
  [tout,yout] = newcorrect('f0', t0, tfinal, h(k-1), y1);
  outputs{k-1} = [tout,yout];
end
for i = 1:(kf-1)
   maxabserror(i) = max(abs(outputs{1,i}(:,2)-f0true(outputs{1,i}(:,1))));
end
 for i = 1:(kf-2)
   consmax(i) = maxabserror(i+1)/maxabserror(i);
end

函数 f0 和 f0true 是：

function yprime = f0(t,y)
yprime = (t.^2 + y.^2)/(2.*t.*y);
end

function y = f0true(t);
y = sqrt(t.*(t+3));
end

但是，当我运行该方法时，问题是尽管为解决方案值提供了非常好的近似值，但无法执行顺序的实验结论（分数变化并且不收敛到$2^{-p}$.

**运行脚本后$kf = 10$，我得到连续的绝对最大误差分数向量：

$$\textbf{consmax} = \big[0.0704, \; 0.0657, \; 0.0639, \; 0.0632, \; 0.0628, \; 0.0627, \; 0.0656, \; 0.2400\big]$$

很容易看出，它的每一项都趋向于这确实是我们想要的，因为是方法的顺序。但是最后一个值是怎么回事？$1/16 = 2^{-4}$$4$

我的错误在哪里？这是可以解释的吗？