The following code display the Predictor-Corrector approach to solve the four dimensional Fractional-Order Delay Differential system,

The below code works fine for some values of $tau$ but it gives an overflow error message for other values of $tau$ for example, when $tau=0.10$ and `n=4500`

(`smax=300`

in the code) I get error message which says that an overflow occurred during computation.

```
a1 = 2.1; b1 = 0.01;
c1 = 2.6;
Subscript(m, 1) = 8.4;
Subscript(m, 2) = 6.4;
Subscript(m, 3) = 2.2;
f(t_, x_, y_, z_, w_) := x*(y - a1) + (Subscript(m, 1)* w) + z;
g(t_, x_, y_, z_, w_) := -(b1*y) + (Subscript(m, 2)* w) - x^2 + 1;
d(t_, x_, y_, z_, w_) := -(c1*z) + (Subscript(m, 3)* w) - x;
q(t_, x_, y_, z_, w_) := -x*y*z;
(Alpha) = 0.93;
h = 0.01;
t(0) = 0;
x0 = 2;
y0 = 0;
z0 = -0.2;
w0 = 0.4;
tau = 0.10; nn = Round(tau/h); smax = 300; n = nn*smax;
Do(x(i) = x0; y(i) = y0; z(i) = z0; w(i) = w0;, {i, -nn, 0});
For(k = 1, k <= n, k++, b(k) = k^(Alpha) - (k - 1)^(Alpha);
a(k) = -(2*k^((Alpha) + 1)) + (k - 1)^((Alpha) + 1) + (k +
1)^((Alpha) + 1););
Do(Do(x1(i) = x(i - nn); y1(i) = y(i - nn);, {i, 0, s*nn});
For(j = (s - 1)*nn + 1, j <= s*nn, j++, t(j) = j*h;
p(j) = (h^(Alpha)*
Sum(b(j - A)*f(A*h, x(A), y1(A), z(A), w(A)), {A, 0, j - 1}))/
Gamma((Alpha) + 1) + x(0);
l(j) = (h^(Alpha)*
Sum(b(j - B)*g(B*h, x1(B), y(B), z(B), w(B)), {B, 0, j - 1}))/
Gamma((Alpha) + 1) + y(0);
r(j) = (h^(Alpha)*
Sum(b(j - C)*d(C*h, x1(C), y(C), z(C), w(C)), {C, 0, j - 1}))/
Gamma((Alpha) + 1) + z(0);
o(j) = (h^(Alpha)*
Sum(b(j - H)*q(H*h, x(H), y(H), z(H), w(H)), {H, 0, j - 1}))/
Gamma((Alpha) + 1) + w(0);
x(j) = (h^(Alpha)*(Sum(
a(j - K)*f(h*K, x(K), y1(K), z(K), w(K)), {K, 1, j - 1}) +
f(h*j, p(j), l(j), r(j),
o(j)) + ((j - 1)^((Alpha) + 1) - (-(Alpha) + j - 1)*
j^(Alpha))*f(0, x(0), y(0), z(0), w(0))))/
Gamma((Alpha) + 2) + x(0);
y(j) = (h^(Alpha)*(Sum(
a(j - F)*g(F*h, x1(F), y(F), z(F), w(F)), {F, 1, j - 1}) +
g(h*j, p(j), l(j), r(j),
o(j)) + ((j - 1)^((Alpha) + 1) - (-(Alpha) + j - 1)*
j^(Alpha))*g(0, x(0), y(0), z(0), w(0))))/
Gamma((Alpha) + 2) + y(0);
z(j) = (h^(Alpha)*(Sum(
a(j - G)*d(G*h, x1(G), y(G), z(G), w(G)), {G, 1, j - 1}) +
d(h*j, p(j), l(j), r(j),
o(j)) + ((j - 1)^((Alpha) + 1) - (-(Alpha) + j - 1)*
j^(Alpha))*d(0, x(0), y(0), z(0), w(0))))/
Gamma((Alpha) + 2) + z(0);
w(j) = (h^(Alpha)*(Sum(
a(j - R)*q(R*h, x1(R), y(R), z(R), w(R)), {R, 1, j - 1}) +
q(h*j, p(j), l(j), r(j),
o(j)) + ((j - 1)^((Alpha) + 1) - (-(Alpha) + j - 1)*
j^(Alpha))*q(0, x(0), y(0), z(0), w(0))))/
Gamma((Alpha) + 2) + w(0););, {s, 1, smax})
lst = Table({t(i), x(i)}, {i, n});
lst1 = Table({x(i), y(i), z(i)}, {i, n});
ListLinePlot(lst, PlotRange -> Full,
Frame -> True, PlotLabel -> "(Alpha)=0.93,(Tau)=0.10",
FrameStyle -> Black,
FrameLabel -> {t, x}, LabelStyle -> Directive(Blue, Bold),
PlotStyle -> Purple)
ListPointPlot3D(lst1, PlotRange -> Full,
AxesLabel -> {x, y, z}, PlotLabel -> "(Alpha)=0.93,(Tau)=0.10",
LabelStyle -> Directive(Black, Bold, FontFamily -> "Helvetica"),
TicksStyle -> Black, PlotStyle -> Purple)
```

Could anyone please solve this issue.