How can I write Matlab code for explicit second derivative two-step Runge Kutta methods? I tried the following code on the stated IVP but works slowly

%IVP;

%y'=f(x,y), y(0) = y_{0};

% y''= f_{y}f= f '(x,y) = g(x,y);

%Example

%f1= -4*y1+3*y2+6; y1(0) = 0;

%f2= -2.4*y1+1.6*y2+3.6; y2(0) = 0;

% g1= -4*f1+3*f2;

% g2= -2.4*f1+1.6*f2;

% Methods:

% c=[1, 0]; Y_{1} ^{[n]} = y(x_{n-1}+c1*h); Y_{1} ^{[n-1]} = y(x_{n-2}+c1*h); Y_{2} ^{[n]} = y(x_{n-1}+c2*h); Y_{2} ^{[n-1]} = y(x_{n-2}+c2*h);

%stage 1;

% Y_{1}=(8/9)+y_{n-1}+(1/9)*y_{n-2}+h*((-52/99)*f(Y_{1}^{[n-1]})+(7/11)*f(Y_{2}^{[n-1]}))+h^2*((707/990)*g(Y_{1}^{[n-1]})+(2/15)*g(Y_{2}^{[n-1]}));

% stage 2;

%Y_{2}=(14/15)+y_{n-1}+(1/15)*y_{n-2}+h*((17/26)*f(Y_{1}^{[n]})+(601/4290)*f(Y_{1}^{[n-1]})+(3/11)*f(Y_{2}^{[n-1]}))...

% +h^2*((-371/85800)*g(Y_{1}^{[n-1]})+(13/20)*g(Y_{2}^{[n-1]}));

%Output method;

%y_{n}= y_{n-1}+h*((1/9)*f(Y_{1}^{[n]})+(1/6)*f(Y_{2}^{[n]})+(13/6)*f(Y_{1}^{[n-1]})-(13/9)*f(Y_{2}^{[n-1]}))...

% +h^2*((1/10)*g(Y_{1}^{[n]})+(1/8)*g(Y_{2}^{[n]})+(7/8)*g(Y_{1}^{[n-1]})+(37/20)*g(Y_{2}^{[n-1]}));

%where n = 2, 3,....

%-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------clear all;

global b

format long e

%starting stepsize

h= 1e-1;

hout=1e-1;

%time

x1= 0.0;

x2= 0.0;

%initial condition

y01= 0.0;

y02= 0.0;

x1out = 0;

x2out = 0;

y1out=0.0;

y2out=0.0;

errout=0;

order

p = 4;

kp = 0.4/( p +1);

kI = 0.3/( p +1);

i=0;

nfe=0;

nrs=0;

tic;

while x1<1;

%first derivative

f01= -4*y01+3*y02+6;

f02= -2.4*y01+1.6*y02+3.6;

%second derivative

g01= -4*f01+3*f02;

g02= -2.4*f01+1.6*f02;

% additional starting value

y11=-3.375*exp(-2*(x1+h))+1.875*exp(-0.4*(x1+h))+1.5;

y12=-2.25*exp(-2*(x2+h))+2.25*exp(-0.4*(x2+h));

f11= -4*y11+3*y12+6;

f12= -2.4*y11+1.6*y12+3.6;

g11= -4*f11+3*f12;

g12= -2.4*f11+1.6*f12;

%stage 1: Y1; c1=0;

Y11=(8/9)+y11+(1/9)*y01+h*((-52/99)*f01+(7/11)*f11)+h^2*((707/990)*g01+(2/15)*g11);

Y12=(8/9)+y12+(1/9)*y02+h*((-52/99)*f02+(7/11)*f12)+h^2*((707/990)*g02+(2/15)*g12);

fY11= -4*Y11+3*Y12+6;

fY12= -2.4*Y11+1.6*Y12+3.6;

gY11= -4*fY11+3*fY12;

gY12= -2.4*fY11+1.6*fY12;

%stage 2: Y2; c2=1;

Y21=(14/15)+y11+(1/15)*y01+h*((17/26)*fY11+(601/4290)*f01+(3/11)*f11)+h^2*((-371/85800)*g01+(13/20)*g11);

Y22=(14/15)+y12+(1/15)*y02++h*((17/26)*fY12+(601/4290)*f02+(3/11)*f12)+h^2*((-371/85800)*g02+(13/20)*g12);

fY21= -4*Y21+3*Y22+6;

fY22= -2.4*Y21+1.6*Y22+3.6;

gY21= -4*fY21+3*fY22;

gY22= -2.4*fY21+1.6*fY22;

%Error estimate

Est=[(h*(49140*f01+50778*fY11+(19048*g01-41769*gY11+10800*gY21)*h))/33306,...

(h*(49140*f02+50778*fY12+(19048*g02-41769*gY12+10800*gY22)*h))/33306];

err=norm(Est);

% Output method

y21= y11+h*((1/9)*fY11+(1/6)*fY21+(13/6)*f01-(13/9)*f11)+h^2*((1/10)*gY11+(1/8)*gY21+(7/8)*g01+(37/20)*g11);

y22= y12+h*((1/9)*fY12+(1/6)*fY22+(13/6)*f02-(13/9)*f12)+h^2*((1/10)*gY12+(1/8)*gY22+(7/8)*g01+(37/20)*g12);

nfe=nfe+4;

%Tolerance

T= 10^(-2);

if (err <= T);

x1=x1+h;

x2=x1;

x1out = [x1out, x1];

x2out = [x2out, x2];

y1out = [y1out, y21];

y2out = [y2out, y22];

errout=[errout, err];

hout = [hout, h];

%exact solution

u1=-3.375*exp(-2*x1out)+1.875*exp(-0.4*x1out)+1.5;

u2=-2.25*exp(-2*x2out)+2.25*exp(-0.4*x2out);

if nrs<=0;

else

nrs=nrs+1;

end

if(err<=1 )

%if(err<=1)

h = h*(T /err) ^(1/( p+1));

else

% New PI step size

h = h*( T/err)^kI *( errold /err)^kp;

end

% Save the error for use in the PI step size controller

errold = err;

i = i+1;

else

% New asymptotic step size

h = 0.9*h*(T /err) ^(1/( p+1));

end

end

plot(x1out, y1out,'ro-','LineWidth',2,'MarkerFaceColor','r','MarkerSize',3)

hold on

plot(x1out,u1,'g-.','LineWidth',2,'MarkerFaceColor','g','MarkerSize',3)

hold on

plot(x2out, y2out,'ko-','LineWidth',2,'MarkerFaceColor','r','MarkerSize',3)

hold on

plot(x2out,u2,'y-.','LineWidth',2,'MarkerFaceColor','g','MarkerSize',3)

hold off

XLABEL('x')

YLABEL('I(x)')

legend('ESDTSRK (y1)','I1(x)','ESDTSRK (y2)', 'I2(x)');

a=[y1out' u1' y2out' u2'];

z1=abs(u1'-y1out');

z2=abs(u2'-y2out');

Z=[z1 z2];

ns=size(x1out)

nrs

nfe

Est

err

toc