can anyone run my code?

please do it and send it for me.

thanks.

restart;

with(plottools): with(LinearAlgebra): with(plots):

ode := `assuming`([diff(Y(X), `$`(X, 2))+2*alpha*(diff(Y(X), X))+beta^2*Y(X) = 0], [alpha >= 0, beta >= 0, alpha+beta > 0]):

F := unapply(solve(subs({X = x, Y(X) = y, diff(Y(X), X) = yp, diff(Y(X), `$`(X, 2)) = yz}, ode), yz), x, y, yp):

Fp := unapply(solve(subs({X = x, Y(X) = y, diff(Y(X), X) = yp, diff(Y(X), `$`(X, 2)) = yz, diff(Y(X), `$`(X, 3)) = yt}, diff(ode, X)), yt), x, y, yp, yz):

Ni := seq(i, i = 0 .. 9), 15, 20:

for i in Ni do

print(ni);

st := time[real]();

f[0, ni] := F(x[0], y[0, ni], yp[0, ni]);

fp[0, ni] := Fp(x[0], y[0, ni], yp[0, ni], f[0, ni]);

f[1, ni] := F(x[1], y[1, ni], yp[1, ni]);

fp[1, ni] := Fp(x[1], y[1, ni], yp[1, ni], f[1, ni]);

y[2, 0] := y[0, ni]+2*h*yp[0, ni]+(6/5)*f[0, ni]*h^2+(4/15)*fp[0, ni]*h^3+(4/5)*f[1, ni]*h^2+(4/15)*fp[1, ni]*h^3;

yp[2, 0] := yp[0, ni]+2*f[0, ni]*h+(2/3)*fp[0, ni]*h^2+(4/3)*fp[1, ni]*h^2;

for j to ni do

f[2, j-1] := F(x[2], y[2, j-1], yp[2, j-1]);

fp[2, j-1] := Fp(x[2], y[2, j-1], yp[2, j-1], f[2, j-1]);

y[2, j] := y[1, ni]+h*yp[1, ni]+(7/20)*f[1, ni]*h^2+(1/20)*fp[1, ni]*h^3+(3/20)*f[2, j-1]*h^2-(1/30)*fp[2, j-1]*h^3;

yp[2, j] := yp[1, ni]+(1/2)*f[1, ni]*h+(1/12)*fp[1, ni]*h^2+(1/2)*f[2, j-1]*h-(1/12)*fp[2, j-1]*h^2;

end do:

Ms := Matrix(4, 4); Ms[1, 3] := 1; Ms[2, 4] := 1;

y[2, ni] := collect(algsubs(h*alpha = H1, expand(algsubs(h*beta = H2, expand(y[2, ni])))), {y[0, ni], y[1, ni], yp[0, ni], yp[1, ni]});

Ms[3, 1] := coeff(y[2, ni], y[0, ni]);

Ms[3, 2] := coeff(y[2, ni], yp[0, ni])/h;

Ms[3, 3] := coeff(y[2, ni], y[1, ni]);

Ms[3, 4] := coeff(y[2, ni], yp[1, ni])/h;

hyp[2, ni] := collect(algsubs(h*alpha = H1, expand(algsubs(h*beta = H2, expand(h*yp[2, ni])))), {y[0, ni], y[1, ni], yp[0, ni], yp[1, ni]});

Ms[4, 1] := coeff(hyp[2, ni], y[0, ni]);

Ms[4, 2] := coeff(hyp[2, ni], yp[0, ni])/h;

Ms[4, 3] := coeff(hyp[2, ni], y[1, ni]);

Ms[4, 4] := coeff(hyp[2, ni], yp[1, ni])/h;

sol := Eigenvalues(Ms);

print(time[real]()-st);

st := time[real]();

SR[ni, 1] := implicitplot(max(seq(abs(sol[ii]), ii = 1 .. numelems(sol))) <= 1, H1 = 0 .. 3, H2 = 0 .. 3, filledregions, gridrefine = 3, axes = Boxed, view = [-2 .. 3, -3 .. 3], labels = [H[1], H[2]], labeldirections = ["horizontal", "vertical"]);

SR[ni, 2] := implicitplot(max(seq(abs(sol[ii]), ii = 1 .. numelems(sol))) <= 1, H1 = -2 .. 3, H2 = -2 .. 3, gridrefine = 3, axes = Boxed, view = [-2 .. 3, -3 .. 3], labels = [H[1], H[2]], labeldirections = ["horizontal", "vertical"]);

print(time[real]()-st);

end do;

for i in Ni do

i;

display({SR[i, 1], SR[i, 2], line([-1, 0], [3, 0], color = red, linestyle = dash), line([0, -3], [0, 3], color = red, linestyle = dash)});

end do;

display({seq(SR[i, 2], i = 0 .. Ni)});