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}
 ]

Ok. Some remarks.

1) I think the equations (1)... are not quite correct, since the "natural" death-rate is described by gg N1. N1 is the total population including the already deceased persons. That means, the already dead die again. According to the paper the number of the dead persons is described by b4[x]. Therefore in my opinion it must be gg ( N1 - b4[x] ). The effect of using gg N1 is that for greater gg fb1 becomes negative, which of course is impossible.

2) I don't know what a fractal fractional plot is. You (perhaps we, but I am afraid I don't have that time) could read the paper line by line to understand what the authors do in fact. And perhaps then it is possible to design some working code.

3) What do you think about writing an email to the authors asking them a) what they did and b) for a code (which could perhaps be transferred to Mathemtica) to obtain those plots in question.

Regards HD

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

Hello Zeeshan,

I would be interested to see that code as well.

Greetings, 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