Or try this

eqs = {
  D[b1[x], x] == -aa b1[x] b2[x] + bb b3[x] - gg N1,
  D[b2[x], x] == aa b1[x] b2[x] - ee b2[x] - dd  b2[x],
  D[b3[x], x] == dd b2[x] - bb b3[x],
  D[b4[x], x] == ee b2[x] + gg N1
  }
ic = {
  b1[0] == 1000,
  b2[0] == 3,
  b3[0] == 0,
  b4[0] == 0
  }

Manipulate[
 vals = {
   aa -> aaa ,
   bb -> bbb ,
   gg -> ggg,
   ee -> eee ,
   dd -> ddd ,
   N1 -> 1000
   };

 sol = Flatten[
   NDSolve[
    Join[eqs, ic] /. vals,
    {b1, b2, b3, b4},
    {x, 0, 30}, MaxSteps -> Infinity
    ]
   ];
 fb1 = b1 /. sol;
 fb2 = b2 /. sol;
 fb3 = b3 /. sol;
 fb4 = b4 /. sol;

 Plot[{fb1[x], fb2[x], fb3[x], fb4[x]}, {x, 0, 20}, 
  PlotRange -> {0, 1000}],

 {{aaa, .004}, .0001, .01},
 {{bbb, .01}, .001, .1},
 {{ggg, .01}, .001, .1},
 {{eee, .02}, .001, .1},
 {{ddd, .6}, .01, 1}
 ]

Hello Zeeshan,

did you get an answer? I spent litterally several hours trying to understand that paper and how they managed to get those plots. But no success. If I had been the referee I would have rejected a paper like this: difficult to understand and moreover no code showing how they programmed the formulas used.

I hope the authors will clarify what they did, and if they haven't answered yet I like to suggest you contact them once more.

Regards HD

But try

vals = {
  aa -> .004 ,
  bb -> .01 ,
  gg -> .0001,
  ee -> .02 ,
  dd -> .6 ,
  N1 -> 1000
  }

this gives something looking similar to the plots you want.

1) Where is your problem? 2) What is your problem?

Long time ago I was told that it is better to avoid indexed variables. I reformulated the 1st part of your code. Does this make sense? b1 and b2 assume very high values. May be this is the reason for the overflow?

   step = 0.02; h = step; tmax = 10; intervals = tmax/h; 
    t[0] = 0.00000001;
    (*Fractional order***************************************)
    \[Rho] = 1;

   (*Fractal order******************************************)
    k = 1;

   (*Parameters*********************************************)\
    \[Alpha] = 0.4; \[Beta] = 0.01; \[Gamma] = 0.01; \[Epsilon] = 0.02; \
    \[Delta] = 0.6;

    (*Ebola equations*****************************************)
    V1[n_] := -\[Alpha]*b1[n]*b2[n] + \[Beta]*b3[n] - \[Gamma]*N1;
    V2[n_] := \[Alpha]*b1[n]*b2[n] + \[Epsilon]*b2[n] - \[Delta]*b2[n];
    V3[n_] := \[Delta]*b2[n] - \[Beta]*b3[n];
    V4[n_] := \[Epsilon]*b2[n] + \[Gamma]*N1;

  (*Total Population****************************************)
    N1 = 1000;

  (*Initial condition***************************************)
    b1[0] = 1000
    b2[0] = 15
    b3[0] = 0
    b4[0] = 0


    (*Fractional \    Euler****************************************)
For[n = 0, n <= 2, n++,
     b1[n + 1] = b1[n] + h^\[Rho]  V1[n]/Gamma[\[Rho] + 1]*(V1[n]);
     b2[n + 1] = b2[n] + h^\[Rho] / Gamma[\[Rho] + 1] V2[n];
     b3[n + 1] = b3[n] + h^\[Rho]/ Gamma[\[Rho] + 1] V3[n];
     b4[n + 1] = b4[n] + h^\[Rho] / Gamma[\[Rho] + 1] V4[n]
     ]







 b1 /@ Range[0, 3]
    b2 /@ Range[0, 3]
    b3 /@ Range[0, 3]
    b4 /@ Range[0, 3]

    Out[16]= {1000, 723402., 3.04409*10^13, 1.80592*10^36}

    Out[17]= {15, 134.826, 780400., 1.90048*10^17}

    Out[18]= {0, 0.18, 1.79788, 9366.6}

    Out[19]= {0, 0.206, 0.45993, 312.82}

Ok.

But there seems to be something severely wrong. Look at the pics: beta1 alias b1 starts at 1000 and decreases to well below 200. But with your (transcribed) code I get

b1[3]
V1[3]

Out[25]= 1.80592*10^36

Out[26]= -1.37285*10^53

Confusing indeed.

I had a look at the formulas around equation (20) and reformulated it. To no avail. There are again very big numbers given back.

And look at eqs (13) to (16). Do they make sense (look at the variables of the functions)?

I am afraid the authors didn't check their paper for typos.

And if I am not wrong their function psi2 always gives 3

(n + 1 - m)^rho (n - m + 2 + rho) - (n - m)^rho (n - m + 2 + 2 rho) // FullSimplify