Wrong results for $2$ stage multistep method $y_{n+2} - y_n = h\left[(1/3)f_{n+2} + (4/3)f_{n+1} + (1/3)f_n\right]$
I need to fix a code to utilise the $2$ stage multistep method :
$$y_{n+2} - y_n = h\left[(1/3)f_{n+2} + (4/3)f_{n+1} + (1/3)f_n\right]$$
Since this is an implicit method, I used a Newton-Raphson approach for the final determination of $y_{n+2}$.
Below, is my attempt at a code implementing the already given rk4 method (runge-kutta 4) for the first approximation steps :
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
This is to be used for showing experimentally that when you divide the step by two, the fraction of the maximum absolute errors per consecutive different steps $h$ is approximately equal to $2^{-p}$ where $p$ is the order of the method. The method is proven to be of order $4$. The exercise whishes the initial starting point to be $k0=2$ thus I created the following script :
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
The function f0 and f0true are :
function yprime = f0(t,y)
yprime = (t.^2 + y.^2)/(2.*t.*y);
end
function y = f0true(t);
y = sqrt(t. 935