I have 4 step Adams - Bashforth method code. It works. But the results are bad. And plot function gives an error

f = @(t, y) cos(2*t + y) + (3/2)*(t - y);
y0 = 1;
a = 0;
b = 1;
N0 = 10;
p = 4; % 4-step Adams-Bashforth
j = 1;
for k = 0:5
    N = N0*2^k;
    NN(j) = N;
    h = (b-a)/N;
    hh(j) = h;
    tn = a:h:b;
    yn = zeros(1,length(tn));
    y(1) = y0;
    % Runge-Kutta for initial values
    for i = 1:p-1
        k1 = h*f(tn(i), y(i));
        k2 = h*f(tn(i) + h/2, y(i) + k1/2);
        k3 = h*f(tn(i) + h/2, y(i) + k2/2);
        k4 = h*f(tn(i) + h, y(i) + k3);
        y(i+1) = y(i) + (k1 + 2*k2 + 2*k3 + k4)/6;
    end
    t = a + (p-1)*h;
    for i = p:length(tn)-1
        % Adams-Bashforth 4-step method
        ynew = fzero(@(ynew) ynew  + (55/24)*y(i) - (59/24)*y(i-1) + ...
            (37/24)*y(i-2) - (9/24)*y(i-3)+(12/24)*h*f(t+h,ynew), y(i));
        y(i+1) = ynew;
        t = t + h;
    end

    plot(tn, y)
    hold on;
    
    yend(j) = y(end);
    j=j+1;
    
end
figure(1);
xlabel('t');
ylabel('y(t)');

[mk, nk] = size(yend);
for k = 1:nk-1
    AB_error(k) = abs(yend(k) - yend(k+1))/(2^p - 1); % Error calculation
end
for k = 1:nk-2
    AB_error_ratio(k) = AB_error(k) / AB_error(k+1); % Ratio calculation
    AB_convergence_speed(k) = log(AB_error(k) / AB_error(k+1)) / log(2); % Convergence order calculation
end

Answer 1

[Torsten]

% Adams-Bashforth 4-step method

ynew = fzero(@(ynew) ynew + (55/24)*y(i) - (59/24)*y(i-1) + ...

(37/24)*y(i-2) - (9/24)*y(i-3)+(12/24)*h*f(t+h,ynew), y(i));

This is not the Adams-Bashford 4-step method. Adams-Bashford 4-step is explicit and given by

y(i+1)= y(i) + h/24*(55*f(tn(i),y(i))-59*f(tn(i-1),y(i-1))+37*f(tn(i-2),y(i-2))-9*f(tn(i-3),y(i-3)))