How can I implement Richardson Extrapolation to the third order Adams-Bashforth Method

So, there is an second order differantial equation of ; which has the initial condition of; . And the analytical solution is;.

the desired local accuracy equal to 10^-4

where the accuracy is defined with:

So we need to expect the condition below;

so far, I have done the analytical solution and Adams-Bashforth 3rd method but could not understand how to implement the richardson method to the Adams method. And I need to get a result of a single plot with two traces: analytical solution and numerical solution satisfying the conditions above.

I know the mathematical formulas of richardson method, but i am not familiar it into matlab. If anyone can help me on this problem, that would be really great.

Thanks all.

These are codes that I have wrote ;

%%% Analytical solution %%%%%

h  = 0.1;
x  = 0:h:5;

c1 = 1;
c2 = (1 - (10*pi)/(1000*(pi)^2 -1 )) / sqrt(10);

y_analytical = c2*sin( sqrt(10).*x ) + c1*cos( sqrt(10).*x )+ sin( 100*pi.*x )/( 10000*pi^2-10 );    

plot (x,y_analytical);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  The code below is a seperate script %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%% ADAM'S BASFORTH METHOD %%%
clear all

h = 0.01;
y = [1;1];
t = 0:h:0.2;

y_analitical = funct_analitical(t);
y(:,2) = y(:,1)+h*funct(t(1),y(:,1));
y(:,3) = y(:,2) + 1/2*h*(3*funct(t(2), y(:,2)) - funct(t(2)-h, y(:,1)));

for i = 3 : length(t)-1
    y_minus2 = [
        spline(t(1:i),y(1,:), t(i)-2*h)
        spline(t(1:i),y(2,:), t(i)-2*h)
    ];
    y_minus1 = [
        spline(t(1:i),y(1,:), t(i)-h)
        spline(t(1:i),y(2,:), t(i)-h)
    ];
    y(:,i+1) = y(:,i) + h*(23/12*funct(t(i), y(:,i)) - 4/3*funct(t(i)-h, y_minus1) +5/12*funct(t(i)-2*h, y_minus2)); % Second order
end

%%% ERROR CALCULATIONS %%% 

avg_error = mean(abs(y(1,:)- y_analitical))
no_of_steps = length(t)
max_error = max (abs(y(1,:)- y_analitical))


figure
plot(t,y_analitical,'black');
hold on
plot(t,y(1,:),'red');
grid on



%%% ANALİTİCAL SOLUTION %%%
function out = funct_analitical(t)

c1 = 1;
c2 = (1 - (10*pi)/(1000*(pi)^2 -1 )) / sqrt(10);

out = c2*sin( sqrt(10).*t ) + c1*cos( sqrt(10).*t ) + sin( 100*pi.*t )/( 10000*pi^2-10 );
end

%%% EQUATİON %%%
function out = funct(t,y)
dy2dt = -10*y(1,1)-sin(100*pi*t);
dy1dt = y(2,1);
out = [dy1dt; dy2dt];
end