Solver for COVID-19 epidemic model with the Caputo fractional derivatives

Question

As it is known in biological system with memory it would be suitable to use fractional derivatives to describe evolution of the system. In a current version of Mathematica 12.1 there is no special solver for integrodifferential equations. Here we show solver with using Haar wavelets for dynamic system (13) presented in a paper M.A. Khan, A. Atangana, Modeling the dynamics of novel coronavirus (2019-nCov) with fractional derivative, Alexandria Eng. J. (2020)

with differential operator replaced with the Caputo definition for fractional derivative as follows $$\frac {d f}{dt}\rightarrow \frac {1}{\Gamma (1-\rho)}\int_0^t{\frac{f'(x)dx}{(t-x)^{\rho}}}$$ The code below allows us to reproduce Figure 7 from the paper linked above. Let define functions

h[x_, k_, m_] := WaveletPsi[HaarWavelet[], m x - k];
h1[x_] := WaveletPhi[HaarWavelet[], x]

Let take $\rho =9/10$, and then we can calculate integrals

Integrate[h[t, k, m], {t, 0, x}, Assumptions -> {k >= 0, m > 0, x > 0}]

Integrate[h1[t], {t, 0, x}, Assumptions -> {x > 0}]

Integrate[h[x, k, m]/(t - x)^(9/10), {x, 0, t}, 
 Assumptions -> {t > 0, k >= 0, m > 0}]

Integrate[h1[x]/(t - x)^(9/10), {x, 0, t}, Assumptions -> {t > 0}]

With these integrals let define functions

p[x_, k_, m_] := Piecewise[{{(1 + k - m*x)/m, k >= 0 && 1/m + (2*k)/m - 2*x < 0 && 
      1/m + k/m - x >= 0 && m > 0}, {(-k + m*x)/m, k >= 0 && 1/m + (2*k)/m - 2*x >= 0 && 
      k/m - x < 0 && 1/m + k/m - x >= 0 && m > 0}}, 0]

p1[x_] := Piecewise[{{1, x > 1}}, x]

pc[t_, k_, m_] := Piecewise[{{10*t^(1/10), k == 0 && 1/m - 2*t >= 0 && m > 0 && t > 0 && 
      1/m + (2*k)/m - 2*t >= 0 && 1/m + k/m - t >= 0}, {(10*(-k + m*t)^(1/10))/m^(1/10), 
     k > 0 && 1/m + (2*k)/m - 2*t >= 0 && k/m - t < 0 && m > 0 && 1/m + k/m - t >= 0}, 
    {(10*((-k + m*t)^(1/10) - 2^(9/10)*(-1 - 2*k + 2*m*t)^(1/10)))/m^(1/10), 
     k > 0 && 1/m + (2*k)/m - 2*t < 0 && 1/m + k/m - t >= 0 && m > 0}, 
    {(10*((-1 - k + m*t)^(1/10) + (-k + m*t)^(1/10) - 2^(9/10)*(-1 - 2*k + 2*m*t)^(1/10)))/
      m^(1/10), k > 0 && 1/m + (2*k)/m - 2*t < 0 && 1/m + k/m - t < 0 && m > 0}, 
    {(5*(2*(m*t)^(1/10) - 2^(9/10)*(-1 + 2*m*t)^(1/10) - 2^(9/10)*(-1 - 2*k + 2*m*t)^(1/10)))/
      m^(1/10), k == 0 && 1/m - 2*t < 0 && 1/m + (2*k)/m - 2*t < 0 && 1/m + k/m - t >= 0 && m > 0}, 
    {(5*(2*(m*t)^(1/10) + 2*(-1 - k + m*t)^(1/10) - 2^(9/10)*(-1 + 2*m*t)^(1/10) - 
        2^(9/10)*(-1 - 2*k + 2*m*t)^(1/10)))/m^(1/10), k == 0 && 1/m - 2*t < 0 && 
      1/m + k/m - t < 0 && m > 0}}, 0]

pc1[t_] := Piecewise[{{-10*((-1 + t)^(1/10) - t^(1/10)), t >= 1}}, 10*t^(1/10)]

Now we have all functions to solve a problem

AbsoluteTiming[ J = 4; M = 2^J; dx = 1/(2*M); 
 Np0 = 8266000; 
  μp (*Natural mortality rate*)= 
  1/(76.79 365); Πp (*Birth rate*)= μp Np0 ; ηp \
(*Contact rate*)= 0.05; ψ (*Transmissibility multiple*) = 
  0.02; ηw (*Disease transmission coeﬃcient*)= 
  0.000001231; θp (*The proportion of asymptomatic \
infection*)= 0.1243; ωp (*Incubation period*)= 
  0.00047876;  ρp (*Incubation period*)= 
  0.005;  τp (*Removal or recovery rate of Ip*)= 
  0.09871;  τap (*Removal or recovery rate of Ap *)= 
  0.854302; ϱp (*Contribution of the virus to M by Ip*)= 
  0.000398; ϖp (*Contribution of the virus to M by Ap*) = 
  0.001; πp(*Removing rate of virus from M*) = 0.01;

 var1 = {Sp1, Ep1, Ip1, Ap1, Rp1, Mp1}; 
 var = {Sp, Ep, Ip, Ap, Rp, Mp}; aco = {aS, aE, aI, aA, aR, aM}; 
 aco1 = {aS1, aE1, aI1, aA1, aR1, aM1}; 
 aco0 = {aS0, aE0, aI0, aA0, aR0, aM0};
 A = 0; xl = Table[A + l dx, {l, 0, 2 M}]; 
 xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, 2 M + 1}];
 Sp1[x_] := 
  Sum[aS[i, j] pc[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aS1 pc1[x]; 
 Sp[x_] := 
  Sum[aS[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aS1 p1[x] + aS0; 
 Ep1[x_] := 
  Sum[aE[i, j] pc[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aE1 pc1[x]; 
 Ep[x_] := 
  Sum[aE[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aE1 p1[x] + aE0; 
 Ip1[x_] := 
  Sum[aI[i, j] pc[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aI1 pc1[x]; 
 Ip[x_] := 
  Sum[aI[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aI1 p1[x] + aI0; 
 Ap1[x_] := 
  Sum[aA[i, j] pc[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aA1 pc1[x]; 
 Ap[x_] := 
  Sum[aA[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aA1 p1[x] + aA0; 
 Rp1[x_] := 
  Sum[aR[i, j] pc[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aR1 pc1[x]; 
 Rp[x_] := 
  Sum[aR[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aR1 p1[x] + aR0; 
 Mp1[x_] := 
  Sum[aM[i, j] pc[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aM1 pc1[x]; 
 Mp[x_] := 
  Sum[aM[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aM1 p1[x] + aM0;

 varM = Join[aco0, aco1, 
   Flatten[Table[{aS[i, j], aE[i, j], aI[i, j], aA[i, j], aR[i, j], 
      aM[i, j]}, {j, 0, J, 1}, {i, 0, 2^j - 1, 1}]]];
 ρ = 9/10; tn = (1/120);
 eq1[t_] := -tn/Gamma[1 - ρ] Sp1[t] + Πp/
    Np0 - μp Sp[t] - ηp Sp[
     t] (Ip[t] + ψ Ap[t])/(Sp[t] + Ep[t] + Ip[t] + Ap[t] + 
       Rp[t]) - Np0 ηw Sp[t] Mp[t]; 
 eq2[t_] := -tn/Gamma[1 - ρ] Ep1[t] + ηp  Sp[
     t] (Ip[t] + ψ Ap[t])/(Sp[t] + Ep[t] + Ip[t] + Ap[t] + 
       Rp[t]) + 
   Np0 ηw Sp[t] Mp[t] - (1 - θp) ωp Ep[
     t] - θp ρp Ep[t] - μp Ep[t];
 eq3[t_] := -tn/Gamma[1 - ρ] Ip1[
     t] + (1 - θp) ωp Ep[t] - (τp + μp) Ip[t]; 
 eq4[t_] := -tn/Gamma[1 - ρ] Ap1[t] + θp ρp Ep[
     t] - (τap + μp) Ap[t]; 
 eq5[t_] := -tn/Gamma[1 - ρ] Rp1[t] + τp Ip[
     t] + τap Ap[t] - μp Rp[t]; 
 eq6[t_] := -tn/Gamma[1 - ρ] Mp1[t] + ϱp Ip[
     t] + ϖp Ap[t] - πp Mp[t];


 eq = Flatten[
   ParallelTable[{eq1[t] == 0, eq2[t] == 0, eq3[t] == 0, eq4[t] == 0, 
     eq5[t] == 0, eq6[t] == 0}, {t, xcol}]];
 Do[icv[i] = {Sp[0] == 8065518/Np0/8 i, Ep[0] == 200000/Np0, 
    Ip[0] == 282/Np0, Ap[0] == 200/Np0, Rp[0] == 0, 
    Mp[0] == 50000/Np0};
  eqM = Join[eq, icv[i]];
  solv[i] = 
   FindRoot[eqM, Table[{varM[[j]], .1}, {j, Length[varM]}], 
    MaxIterations -> 1000];
  lstSv[i] = 
   Table[{x 120 , Np0 Evaluate[Sp[x] /. solv[i]]}, {x, 0, 1, .01}]; 
  lstEv[i] = 
   Table[{x 120, Np0 Evaluate[Ep[x] /. solv[i]]}, {x, 0, 1, .01}]; 
  lstIv[i] = 
   Table[{x 120, Np0 Evaluate[Ip[x] /. solv[i]]}, {x, 0, 1, .01}]; 
  lstAv[i] = 
   Table[{x 120, Np0 Evaluate[Ap[x] /. solv[i]]}, {x, 0, 1, .01}]; 
  lstRv[i] = 
   Table[{x 120, Np0 Evaluate[Rp[x] /. solv[i]]}, {x, 0, 1, .01}]; 
  lstMv[i] = 
   Table[{x 120, Np0 Evaluate[Mp[x] /. solv[i]]}, {x, 0, 
     1, .01}];, {i, 1, 8}]]

Finally we visualize solution

{ListLinePlot[Table[lstSv[i], {i, 1, 8}], Frame -> True, 
     FrameLabel -> {"t, days", "\!\(\*SubscriptBox[\(S\), \(p\)]\)"}, 
  PlotRange -> All], 
   ListLinePlot[Table[lstEv[i], {i, 1, 8}], Frame -> True, 
     FrameLabel -> {"t, days", "\!\(\*SubscriptBox[\(E\), \(p\)]\)"}, 
  PlotRange -> All], 
   ListLinePlot[Table[lstIv[i], {i, 1, 8}], Frame -> True, 
     FrameLabel -> {"t, days", "\!\(\*SubscriptBox[\(I\), \(p\)]\)"}, 
  PlotRange -> All], 
   ListLinePlot[Table[lstAv[i], {i, 1, 8}], Frame -> True, 
     FrameLabel -> {"t, days", "\!\(\*SubscriptBox[\(A\), \(p\)]\)"}, 
  PlotRange -> All], 
   ListLinePlot[Table[lstRv[i], {i, 1, 8}], Frame -> True, 
     FrameLabel -> {"t, days", "\!\(\*SubscriptBox[\(R\), \(p\)]\)"}, 
  PlotRange -> All], 
   ListLinePlot[Table[lstMv[i], {i, 1, 8}], Frame -> True, 
  FrameLabel -> {"t, days", "M"}, 
     PlotRange -> All, PlotLegends -> Automatic]}

The question is about how to add $\rho$ in this code as parameter?

Update 1. The straight forward solution of this problem is simply to include $\rho$ in the pc, pc1functions definitions as follows (here $\rho$ is replaced by q) :

pc[t_, k_, m_, q_] := 
 Piecewise[{{-(t^(1 - q)/(-1 + q)), k == 0 && 1/m - 2*t >= 0 && 
          m > 0 && t > 0 && 1/m - t >= 0}, 
      {-((m^(-1 + q)*(1/(-k + m*t))^(-1 + q))/(-1 + q)), 
        k > 0 && 1/m + (2*k)/m - 2*t > 0 && k/m - t < 0 && m > 0 && 
          1/m + k/m - t > 0}, 
      {(-t^q + 2*m*t^(1 + q) - m*t*(-(1/(2*m)) + t)^q)/
          (t^q*(-(1/(2*m)) + t)^q*(m*(-1 + q))), 
        k == 0 && m > 0 && 1/m - 2*t < 0 && 1/m - t >= 0}, 
      {(1/(-1 + q))*((2^(-1 + q)*m^(-1 + 2*q)*(-(-(k/m) + t)^q - 
                   2*k*(-(k/m) + t)^q + 2*m*t*(-(k/m) + t)^q + 
                   2*k*(-((1/2 + k)/m) + t)^q - 
           2*m*t*(-((1/2 + k)/m) + t)^
                       q))/((1 + 2*k - 2*m*t)*(k - m*t))^q), 
        k > 0 && 1/m + (2*k)/m - 2*t == 0 && m > 0 && 
          1/m + k/m - t > 0}, 
      {-((1/(-1 + q))*((2^(-1 + q)*m^(-1 + 2*q)*
                   (-2*(-((1/2 + k)/m) + t)^
               q*((1 + 2*k - 2*m*t)*(k - m*t))^
                          q - 2*k*(-((1/2 + k)/m) + t)^q*
                        ((1 + 2*k - 2*m*t)*(k - m*t))^q + 
                      2*m*t*(-((1/2 + k)/m) + t)^q*((1 + 2*k - 2*m*t)*
                             (k - m*t))^q + (-((1 + k)/m) + t)^q*
                        ((1 + 2*k - 2*m*t)*(k - m*t))^q + 

             2*k*(-((1 + k)/m) + t)^q*((1 + 2*k - 2*m*t)*(k - m*t))^
                          q - 2*m*t*(-((1 + k)/m) + t)^q*
                        ((1 + 2*k - 2*m*t)*(k - m*t))^
               q + (-(k/m) + t)^q*
                        ((1 + 2*k - 2*m*t)*(1 + k - m*t))^q + 

             2*k*(-(k/m) + t)^q*((1 + 2*k - 2*m*t)*(1 + k - m*t))^q - 

             2*m*t*(-(k/m) + t)^q*((1 + 2*k - 2*m*t)*(1 + k - m*t))^
                          q - 2*k*(-((1/2 + k)/m) + t)^q*
                        ((1 + 2*k - 2*m*t)*(1 + k - m*t))^q + 
                      2*m*t*(-((1/2 + k)/m) + t)^q*((1 + 2*k - 2*m*t)*
                             (1 + k - m*t))^
               q))/(((1 + 2*k - 2*m*t)*(k - m*t))^q*
                   ((1 + 2*k - 2*m*t)*(1 + k - m*t))^q))), 
        k > 0 && m > 0 && 1/m + (2*k)/m - 2*t <= 0 && 
          1/m + k/m - t <= 0}, 
      {-((1/(2*m*(-1 + q)))*((2^q*m^(2*q)*t^q*(-(1/m) + t)^q*
                     (-(1/(2*m)) + t)^q - 
           2^(1 + q)*m^(1 + 2*q)*t^(1 + q)*
                     (-(1/m) + t)^q*(-(1/(2*m)) + t)^q - 
           2^(1 + q)*m^(2*q)*
                     t^q*(-(1/(2*m)) + t)^(2*q) + 
           2^(1 + q)*m^(1 + 2*q)*
                     t^(1 + q)*(-(1/(2*m)) + t)^(2*q) + 
                   t^q*((-1 + m*t)*(-1 + 2*m*t))^q - 2*m*t^(1 + q)*
                     ((-1 + m*t)*(-1 + 2*m*t))^q + 
           2*m*t*(-(1/(2*m)) + t)^q*
                     ((-1 + m*t)*(-1 + 2*m*t))^q)/(t^
            q*(-(1/(2*m)) + t)^q*
                   ((-1 + m*t)*(-1 + 2*m*t))^q))), 
        k == 0 && 1/m - 2*t < 0 && 1/m - t < 0 && m > 0}, 
      {(1/(-1 + q))*((2^(-1 + q)*m^(-1 + q)*((-m^q)*(-(k/m) + t)^q - 
                   2*k*m^q*(-(k/m) + t)^q + 
           2*m^(1 + q)*t*(-(k/m) + t)^q + 
                   2*k*m^q*(-((1/2 + k)/m) + t)^q - 2*m^(1 + q)*t*
                     (-((1/2 + k)/m) + t)^
             q - ((1 + 2*k - 2*m*t)*(k - m*t))^q*
                     (1/(-1 - 2*k + 2*m*t))^q - 
                   2*k*((1 + 2*k - 2*m*t)*(k - m*t))^q*
                     (1/(-1 - 2*k + 2*m*t))^q + 
                   2*m*t*((1 + 2*k - 2*m*t)*(k - m*t))^q*
                     (1/(-1 - 2*k + 2*m*t))^q))/((1 + 2*k - 
            2*m*t)*(k - m*t))^
               q), 1/m + (2*k)/m - 2*t < 0 && k > 0 && m > 0 && 
          1/m + k/m - t > 0}}, 0]

pc1[t_, q_] := Piecewise[{{-(t^(1 - q)/(-1 + q)), t <= 1}}, 
    -(((-1 + t)^q*t + t^q - t^(1 + q))/((-1 + t)^q*t^q*(-1 + q)))]

With these functions we can calculate Figure 6 from the paper above with the next piece of code

AbsoluteTiming[J = 4; M = 2^J; dx = 1/(2*M); 
 Np0 = 8266000; 
  \[Mu]p (*Natural mortality rate*)= 
  1/(76.79 365); \[CapitalPi]p (*Birth rate*)= \[Mu]p Np0 ; \[Eta]p \
(*Contact rate*)= 0.05; \[Psi] (*Transmissibility multiple*) = 
  0.02; \[Eta]w (*Disease transmission coeﬃcient*)= 
  0.000001231; \[Theta]p (*The proportion of asymptomatic \
infection*)= 0.1243; \[Omega]p (*Incubation period*)= 
  0.00047876;  \[Rho]p (*Incubation period*)= 
  0.005;  \[Tau]p (*Removal or recovery rate of Ip*)= 
  0.09871;  \[Tau]ap (*Removal or recovery rate of Ap *)= 
  0.854302; \[CurlyRho]p (*Contribution of the virus to M by Ip*)= 
  0.000398; \[CurlyPi]p (*Contribution of the virus to M by Ap*) = 
  0.001; \[Pi]p(*Removing rate of virus from M*) = 0.01;

 var1 = {Sp1, Ep1, Ip1, Ap1, Rp1, Mp1}; 
 var = {Sp, Ep, Ip, Ap, Rp, Mp}; aco = {aS, aE, aI, aA, aR, aM}; 
 aco1 = {aS1, aE1, aI1, aA1, aR1, aM1}; 
 aco0 = {aS0, aE0, aI0, aA0, aR0, aM0};
 A = 0; xl = Table[A + l dx, {l, 0, 2 M}]; 
 xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, 2 M + 1}];
 Sp1[x_, q_] := 
  Sum[aS[i, j] pc[x, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aS1 pc1[x, q]; 
 Sp[x_] := 
  Sum[aS[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aS1 p1[x] + aS0; 
 Ep1[x_, q_] := 
  Sum[aE[i, j] pc[x, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aE1 pc1[x, q]; 
 Ep[x_] := 
  Sum[aE[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aE1 p1[x] + aE0; 
 Ip1[x_, q_] := 
  Sum[aI[i, j] pc[x, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aI1 pc1[x, q]; 
 Ip[x_] := 
  Sum[aI[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aI1 p1[x] + aI0; 
 Ap1[x_, q_] := 
  Sum[aA[i, j] pc[x, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aA1 pc1[x, q]; 
 Ap[x_] := 
  Sum[aA[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aA1 p1[x] + aA0; 
 Rp1[x_, q_] := 
  Sum[aR[i, j] pc[x, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aR1 pc1[x, q]; 
 Rp[x_] := 
  Sum[aR[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aR1 p1[x] + aR0; 
 Mp1[x_, q_] := 
  Sum[aM[i, j] pc[x, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aM1 pc1[x, q]; 
 Mp[x_] := 
  Sum[aM[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aM1 p1[x] + aM0;

 varM = Join[aco0, aco1, 
   Flatten[Table[{aS[i, j], aE[i, j], aI[i, j], aA[i, j], aR[i, j], 
      aM[i, j]}, {j, 0, J, 1}, {i, 0, 2^j - 1, 1}]]];
 tn[q_] := (1/120)^q;
 eq1[t_, q_] := -tn[q]/Gamma[1 - q] Sp1[t, q] + \[CapitalPi]p/
    Np0 - \[Mu]p Sp[t] - \[Eta]p Sp[
     t] (Ip[t] + \[Psi] Ap[t])/(Sp[t] + Ep[t] + Ip[t] + Ap[t] + 
       Rp[t]) - Np0 \[Eta]w Sp[t] Mp[t]; 
 eq2[t_, q_] := -tn[q]/Gamma[1 - q] Ep1[t, q] + \[Eta]p  Sp[
     t] (Ip[t] + \[Psi] Ap[t])/(Sp[t] + Ep[t] + Ip[t] + Ap[t] + 
       Rp[t]) + 
   Np0 \[Eta]w Sp[t] Mp[t] - (1 - \[Theta]p) \[Omega]p Ep[
     t] - \[Theta]p \[Rho]p Ep[t] - \[Mu]p Ep[t];
 eq3[t_, q_] := -tn[q]/Gamma[1 - q] Ip1[t, 
     q] + (1 - \[Theta]p) \[Omega]p Ep[t] - (\[Tau]p + \[Mu]p) Ip[t]; 
 eq4[t_, q_] := -tn[q]/Gamma[1 - q] Ap1[t, q] + \[Theta]p \[Rho]p Ep[
     t] - (\[Tau]ap + \[Mu]p) Ap[t]; 
 eq5[t_, q_] := -tn[q]/Gamma[1 - q] Rp1[t, q] + \[Tau]p Ip[
     t] + \[Tau]ap Ap[t] - \[Mu]p Rp[t]; 
 eq6[t_, q_] := -tn[q]/Gamma[1 - q] Mp1[t, q] + \[CurlyRho]p Ip[
     t] + \[CurlyPi]p Ap[t] - \[Pi]p Mp[t];


 eq[q_] := 
  Flatten[ParallelTable[{eq1[t, q] == 0, eq2[t, q] == 0, 
     eq3[t, q] == 0, eq4[t, q] == 0, eq5[t, q] == 0, 
     eq6[t, q] == 0}, {t, xcol}]];
 Do[icv[i] = {Sp[0] == 8065518/Np0, Ep[0] == 200000/Np0, 
    Ip[0] == 282/Np0, Ap[0] == 200/Np0, Rp[0] == 0, 
    Mp[0] == 50000/Np0};
  eqM[i] = Join[eq[i], icv[i]];
  solv[i] = 
   FindRoot[eqM[i], Table[{varM[[j]], .1}, {j, Length[varM]}], 
    MaxIterations -> 1000];
  lstSv[i] = 
   Table[{x 120 , Np0 Evaluate[Sp[x] /. solv[i]]}, {x, 0, 1, .01}]; 
  lstEv[i] = 
   Table[{x 120, Np0 Evaluate[Ep[x] /. solv[i]]}, {x, 0, 1, .01}]; 
  lstIv[i] = 
   Table[{x 120, Np0 Evaluate[Ip[x] /. solv[i]]}, {x, 0, 1, .01}]; 
  lstAv[i] = 
   Table[{x 120, Np0 Evaluate[Ap[x] /. solv[i]]}, {x, 0, 1, .01}]; 
  lstRv[i] = 
   Table[{x 120, Np0 Evaluate[Rp[x] /. solv[i]]}, {x, 0, 1, .01}]; 
  lstMv[i] = 
   Table[{x 120, Np0 Evaluate[Mp[x] /. solv[i]]}, {x, 0, 
     1, .01}];, {i, {99/100, 9/10, 8/10, 7/10, 6/10}}];]

We can check that it is run 4-5 times longer than code with a fixed $\rho$. Visualization:

{ListLinePlot[Table[lstSv[i], {i, {99/100, 9/10, 8/10, 7/10, 6/10}}], 
  Frame -> True, 
     FrameLabel -> {"t, days", "\!\(\*SubscriptBox[\(S\), \(p\)]\)"}, 
  PlotRange -> All], 
   ListLinePlot[
  Table[lstEv[i], {i, {99/100, 9/10, 8/10, 7/10, 6/10}}], 
  Frame -> True, 
     FrameLabel -> {"t, days", "\!\(\*SubscriptBox[\(E\), \(p\)]\)"}, 
  PlotRange -> All], 
   ListLinePlot[
  Table[lstIv[i], {i, {99/100, 9/10, 8/10, 7/10, 6/10}}], 
  Frame -> True, 
     FrameLabel -> {"t, days", "\!\(\*SubscriptBox[\(I\), \(p\)]\)"}, 
  PlotRange -> All], 
   ListLinePlot[
  Table[lstAv[i], {i, {99/100, 9/10, 8/10, 7/10, 6/10}}], 
  Frame -> True, 
     FrameLabel -> {"t, days", "\!\(\*SubscriptBox[\(A\), \(p\)]\)"}, 
  PlotRange -> All], 
   ListLinePlot[
  Table[lstRv[i], {i, {99/100, 9/10, 8/10, 7/10, 6/10}}], 
  Frame -> True, 
     FrameLabel -> {"t, days", "\!\(\*SubscriptBox[\(R\), \(p\)]\)"}, 
  PlotRange -> All], 
   ListLinePlot[
  Table[lstMv[i], {i, {99/100, 9/10, 8/10, 7/10, 6/10}}], 
  Frame -> True, FrameLabel -> {"t, days", "M"}, 
     PlotRange -> All, PlotLegends -> Automatic]}

Update 2. We can reduce time by 3-4 times simply replace where it is possible function definition f[x_,...]:=... with f=Compile[{{x,_Real},{...}},...]. So in the last code we have to replace first part as follows

h = Compile[{{x, _Real}, {k, _Integer}, {m, _Integer}}, 
   WaveletPsi[HaarWavelet[], m x - k]];

p = Compile[{{x, _Real}, {k, _Integer}, {m, _Integer}}, 
   Piecewise[{{(1 + k - m*x)/m, k >= 0 && 1/m + (2*k)/m - 2*x < 0 && 
             1/m + k/m - x >= 0 && m > 0}, {(-k + m*x)/m, 
      k >= 0 && 1/m + (2*k)/m - 2*x >= 0 && 
             k/m - x < 0 && 1/m + k/m - x >= 0 && m > 0}}, 0]];
h1 = Compile[{{x, _Real}}, WaveletPhi[HaarWavelet[], x]];

p1 = Compile[{{x, _Real}}, Piecewise[{{1, x > 1}}, x]];

pc = Compile[{{t, _Real}, {k, _Integer}, {m, _Integer}, {q, _Real}}, 
   Piecewise[{{-(t^(1 - q)/(-1 + q)), k == 0 && 1/m - 2*t >= 0 && 
            m > 0 && t > 0 && 1/m - t >= 0}, 
        {-((m^(-1 + q)*(1/(-k + m*t))^(-1 + q))/(-1 + q)), 
          k > 0 && 1/m + (2*k)/m - 2*t > 0 && k/m - t < 0 && m > 0 && 
            1/m + k/m - t > 0}, 
        {(-t^q + 2*m*t^(1 + q) - m*t*(-(1/(2*m)) + t)^q)/
            (t^q*(-(1/(2*m)) + t)^q*(m*(-1 + q))), 
          k == 0 && m > 0 && 1/m - 2*t < 0 && 1/m - t >= 0}, 
        {(1/(-1 + q))*((2^(-1 + q)*m^(-1 + 2*q)*(-(-(k/m) + t)^q - 
                     2*k*(-(k/m) + t)^q + 2*m*t*(-(k/m) + t)^q + 
                     2*k*(-((1/2 + k)/m) + t)^q - 
             2*m*t*(-((1/2 + k)/m) + t)^
                         q))/((1 + 2*k - 2*m*t)*(k - m*t))^q), 
          k > 0 && 1/m + (2*k)/m - 2*t == 0 && m > 0 && 
            1/m + k/m - t > 0}, 
        {-((1/(-1 + q))*((2^(-1 + q)*m^(-1 + 2*q)*
                     (-2*(-((1/2 + k)/m) + t)^
                 q*((1 + 2*k - 2*m*t)*(k - m*t))^
                            q - 2*k*(-((1/2 + k)/m) + t)^q*
                          ((1 + 2*k - 2*m*t)*(k - m*t))^q + 

               2*m*t*(-((1/2 + k)/m) + t)^q*((1 + 2*k - 2*m*t)*
                               (k - m*t))^q + (-((1 + k)/m) + t)^q*
                          ((1 + 2*k - 2*m*t)*(k - m*t))^q + 

               2*k*(-((1 + k)/m) + t)^q*((1 + 2*k - 2*m*t)*(k - m*t))^
                            q - 2*m*t*(-((1 + k)/m) + t)^q*
                          ((1 + 2*k - 2*m*t)*(k - m*t))^
                 q + (-(k/m) + t)^q*
                          ((1 + 2*k - 2*m*t)*(1 + k - m*t))^q + 

               2*k*(-(k/m) + t)^q*((1 + 2*k - 2*m*t)*(1 + k - m*t))^
                 q - 

               2*m*t*(-(k/m) + t)^q*((1 + 2*k - 2*m*t)*(1 + k - m*t))^
                            q - 2*k*(-((1/2 + k)/m) + t)^q*
                          ((1 + 2*k - 2*m*t)*(1 + k - m*t))^q + 

               2*m*t*(-((1/2 + k)/m) + t)^q*((1 + 2*k - 2*m*t)*
                               (1 + k - m*t))^
                 q))/(((1 + 2*k - 2*m*t)*(k - m*t))^q*
                     ((1 + 2*k - 2*m*t)*(1 + k - m*t))^q))), 
          k > 0 && m > 0 && 1/m + (2*k)/m - 2*t <= 0 && 
            1/m + k/m - t <= 0}, 
        {-((1/(2*m*(-1 + q)))*((2^q*m^(2*q)*t^q*(-(1/m) + t)^q*
                       (-(1/(2*m)) + t)^q - 
             2^(1 + q)*m^(1 + 2*q)*t^(1 + q)*
                       (-(1/m) + t)^q*(-(1/(2*m)) + t)^q - 
             2^(1 + q)*m^(2*q)*
                       t^q*(-(1/(2*m)) + t)^(2*q) + 
             2^(1 + q)*m^(1 + 2*q)*
                       t^(1 + q)*(-(1/(2*m)) + t)^(2*q) + 
                     t^q*((-1 + m*t)*(-1 + 2*m*t))^q - 2*m*t^(1 + q)*
                       ((-1 + m*t)*(-1 + 2*m*t))^q + 
             2*m*t*(-(1/(2*m)) + t)^q*
                       ((-1 + m*t)*(-1 + 2*m*t))^q)/(t^
              q*(-(1/(2*m)) + t)^q*
                     ((-1 + m*t)*(-1 + 2*m*t))^q))), 
          k == 0 && 1/m - 2*t < 0 && 1/m - t < 0 && m > 0}, 
        {(1/(-1 + q))*((2^(-1 + q)*m^(-1 + q)*((-m^q)*(-(k/m) + t)^q - 
                     2*k*m^q*(-(k/m) + t)^q + 
             2*m^(1 + q)*t*(-(k/m) + t)^q + 
                     2*k*m^q*(-((1/2 + k)/m) + t)^q - 2*m^(1 + q)*t*
                       (-((1/2 + k)/m) + t)^
               q - ((1 + 2*k - 2*m*t)*(k - m*t))^q*
                       (1/(-1 - 2*k + 2*m*t))^q - 
                     2*k*((1 + 2*k - 2*m*t)*(k - m*t))^q*
                       (1/(-1 - 2*k + 2*m*t))^q + 
                     2*m*t*((1 + 2*k - 2*m*t)*(k - m*t))^q*
                       (1/(-1 - 2*k + 2*m*t))^q))/((1 + 2*k - 
              2*m*t)*(k - m*t))^
                 q), 1/m + (2*k)/m - 2*t < 0 && k > 0 && m > 0 && 
            1/m + k/m - t > 0}}, 0]];

pc1 = Compile[{{t, _Real}, {q, _Real}}, 
  Piecewise[{{-(t^(1 - q)/(-1 + q)), t <= 1}}, 
     -(((-1 + t)^q*t + t^q - t^(1 + q))/((-1 + t)^q*
        t^q*(-1 + q)))]]; tn = Compile[{{q, _Real}}, (1/120)^q];

"it would be suitable to use fractional derivatives to describe evolution of the system" - just out of curiosity, do you have a good reference that demonstrates this? — Chris K, Commented May 10, 2020 at 19:12
@Szabolcs This code is quite good for particular $\rho$ (it takes few second to calculate solution). To make it faster I have introduced functions like pc[t,k,m], pc1[t,k,m] calculated with a given $\rho$. But when I try to use it in general form as a functions of $\rho$, it takes hundred times longer. — Alex Trounev, Commented May 10, 2020 at 20:09
@ChrisK I have few papers about application of fractional derivatives to biological systems including this one arxiv.org/abs/2005.01820 — Alex Trounev, Commented May 10, 2020 at 20:28
Note that there is now ResourceFunction["FractionalIntegrate"] in the Wolfram Function Repository. — Daniel Lichtblau, Commented May 11, 2020 at 15:57
@Daniel, I think Alex would use ResourceFunction["FractionalD"] instead for this, since he wants the Caputo derivative. — J. M.'s missing motivation, Commented May 11, 2020 at 16:54

Rolf Mertig · Accepted Answer · 2020-05-18 13:42:38Z

WaveletPsi is not compilable, so I changed some of your compile definitions. Changing the parallel evaluation a little bit your code runs now in about a second on a simple 4-core machine. Using the experimental FunctionCompile does not help at all, since both the Piecewise and the HaarWavelet function generate compile errors. Looks like those functions are not among the "approximately 2000 functions which cover 31 functionality areas" [Compiler paper].

h = Function[{x, k, m}, WaveletPsi[HaarWavelet[], m x - k]]; 

p = Compile[{{x, _Real}, {k, _Integer}, {m, _Integer}}, 
   Piecewise[{{(1 + k - m*x)/m, 
      k >= 0 && 1/m + (2*k)/m - 2*x < 0 && 1/m + k/m - x >= 0 && 
       m > 0}, {(-k + m*x)/m, 
      k >= 0 && 1/m + (2*k)/m - 2*x >= 0 && k/m - x < 0 && 
       1/m + k/m - x >= 0 && m > 0}}, 0]];

h1 = Function[{x}, WaveletPhi[HaarWavelet[], x]];

p1 = Function[x, Piecewise[{{1, x > 1}}, x]];

pc = Compile[{{t, _Real}, {k, _Integer}, {m, _Integer}, {q, _Real}}, 
   Piecewise[{{-(t^(1 - q)/(-1 + q)), 
      k == 0 && 1/m - 2*t >= 0 && m > 0 && t > 0 && 
       1/m - t >= 
        0}, {-((m^(-1 + q)*(1/(-k + m*t))^(-1 + q))/(-1 + q)), 
      k > 0 && 1/m + (2*k)/m - 2*t > 0 && k/m - t < 0 && m > 0 && 
       1/m + k/m - t > 
        0}, {(-t^q + 2*m*t^(1 + q) - 
         m*t*(-(1/(2*m)) + t)^q)/(t^q*(-(1/(2*m)) + t)^
          q*(m*(-1 + q))), 
      k == 0 && m > 0 && 1/m - 2*t < 0 && 
       1/m - t >= 
        0}, {(1/(-1 + q))*((2^(-1 + q)*
           m^(-1 + 2*q)*(-(-(k/m) + t)^q - 2*k*(-(k/m) + t)^q + 
             2*m*t*(-(k/m) + t)^q + 2*k*(-((1/2 + k)/m) + t)^q - 
             2*m*t*(-((1/2 + k)/m) + t)^q))/((1 + 2*k - 2*m*t)*(k - 
              m*t))^q), 
      k > 0 && 1/m + (2*k)/m - 2*t == 0 && m > 0 && 
       1/m + k/m - t > 
        0}, {-((1/(-1 + q))*((2^(-1 + q)*
             m^(-1 + 
                2*q)*(-2*(-((1/2 + k)/m) + t)^
                 q*((1 + 2*k - 2*m*t)*(k - m*t))^q - 
               2*k*(-((1/2 + k)/m) + t)^
                 q*((1 + 2*k - 2*m*t)*(k - m*t))^q + 
               2*m*t*(-((1/2 + k)/m) + t)^
                 q*((1 + 2*k - 2*m*t)*(k - m*t))^
                 q + (-((1 + k)/m) + t)^
                 q*((1 + 2*k - 2*m*t)*(k - m*t))^q + 
               2*k*(-((1 + k)/m) + t)^q*((1 + 2*k - 2*m*t)*(k - m*t))^
                 q - 2*m*
                t*(-((1 + k)/m) + t)^q*((1 + 2*k - 2*m*t)*(k - m*t))^
                 q + (-(k/m) + t)^q*((1 + 2*k - 2*m*t)*(1 + k - m*t))^
                 q + 2*k*(-(k/m) + t)^
                 q*((1 + 2*k - 2*m*t)*(1 + k - m*t))^q - 
               2*m*t*(-(k/m) + t)^q*((1 + 2*k - 2*m*t)*(1 + k - m*t))^
                 q - 2*k*(-((1/2 + k)/m) + t)^
                 q*((1 + 2*k - 2*m*t)*(1 + k - m*t))^q + 
               2*m*t*(-((1/2 + k)/m) + t)^
                 q*((1 + 2*k - 2*m*t)*(1 + k - m*t))^
                 q))/(((1 + 2*k - 2*m*t)*(k - m*t))^
              q*((1 + 2*k - 2*m*t)*(1 + k - m*t))^q))), 
      k > 0 && m > 0 && 1/m + (2*k)/m - 2*t <= 0 && 
       1/m + k/m - t <= 
        0}, {-((1/(2*
             m*(-1 + q)))*((2^q*m^(2*q)*
              t^q*(-(1/m) + t)^q*(-(1/(2*m)) + t)^q - 
             2^(1 + q)*m^(1 + 2*q)*
              t^(1 + q)*(-(1/m) + t)^q*(-(1/(2*m)) + t)^q - 
             2^(1 + q)*m^(2*q)*t^q*(-(1/(2*m)) + t)^(2*q) + 
             2^(1 + q)*m^(1 + 2*q)*t^(1 + q)*(-(1/(2*m)) + t)^(2*q) + 
             t^q*((-1 + m*t)*(-1 + 2*m*t))^q - 
             2*m*t^(1 + q)*((-1 + m*t)*(-1 + 2*m*t))^q + 
             2*m*t*(-(1/(2*m)) + t)^q*((-1 + m*t)*(-1 + 2*m*t))^q)/(t^
              q*(-(1/(2*m)) + t)^q*((-1 + m*t)*(-1 + 2*m*t))^q))), 
      k == 0 && 1/m - 2*t < 0 && 1/m - t < 0 && 
       m > 0}, {(1/(-1 + q))*((2^(-1 + q)*
           m^(-1 + q)*((-m^q)*(-(k/m) + t)^q - 
             2*k*m^q*(-(k/m) + t)^q + 2*m^(1 + q)*t*(-(k/m) + t)^q + 
             2*k*m^q*(-((1/2 + k)/m) + t)^q - 
             2*m^(1 + q)*
              t*(-((1/2 + k)/m) + t)^
               q - ((1 + 2*k - 2*m*t)*(k - m*t))^
               q*(1/(-1 - 2*k + 2*m*t))^q - 
             2*k*((1 + 2*k - 2*m*t)*(k - m*t))^
               q*(1/(-1 - 2*k + 2*m*t))^q + 
             2*m*t*((1 + 2*k - 2*m*t)*(k - m*t))^
               q*(1/(-1 - 2*k + 2*m*t))^q))/((1 + 2*k - 2*m*t)*(k - 
              m*t))^q), 
      1/m + (2*k)/m - 2*t < 0 && k > 0 && m > 0 && 
       1/m + k/m - t > 0}}, 0]];

pc1 = Compile[{{t, _Real}, {q, _Real}}, 
   Piecewise[{{-(t^(1 - q)/(-1 + q)), 
      t <= 1}}, -(((-1 + t)^q*t + t^q - t^(1 + q))/((-1 + t)^q*
         t^q*(-1 + q)))]];

tn = Function[{q}, (1/120)^q];


(*Now we have all functions to solve a problem with the given \
parameres*)

LaunchKernels[] (* Launch parallel kernels before doing \
AbsoluteTiming *)

AbsoluteTiming[
 J = 4; M = 2^J; dx = 1/(2*M);
 Np0 = 8266000;
 \[Mu]p (*Natural mortality rate*)= 
  1/(76.79 365); \[CapitalPi]p (*Birth rate*)= \[Mu]p Np0; \[Eta]p \
(*Contact rate*)= 0.05; \[Psi] (*Transmissibility multiple*)= 
  0.02; \[Eta]w (*Disease transmission coeﬃcient*)= 
  0.000001231; \[Theta]p (*The proportion of asymptomatic infection*)=
   0.1243; \[Omega]p (*Incubation period*)= 
  0.00047876; \[Rho]p (*Incubation period*)= 
  0.005; \[Tau]p (*Removal or recovery rate of Ip*)= 
  0.09871; \[Tau]ap (*Removal or recovery rate of Ap*)= 
  0.854302; \[CurlyRho]p (*Contribution of the virus to M by Ip*)= 
  0.000398; \[CurlyPi]p (*Contribution of the virus to M by Ap*)= 
  0.001; \[Pi]p(*Removing rate of virus from M*)= 0.01;
 var1 = {Sp1, Ep1, Ip1, Ap1, Rp1, Mp1};
 var = {Sp, Ep, Ip, Ap, Rp, Mp}; aco = {aS, aE, aI, aA, aR, aM};
 aco1 = {aS1, aE1, aI1, aA1, aR1, aM1};
 aco0 = {aS0, aE0, aI0, aA0, aR0, aM0};
 A = 0; xl = Table[A + l dx, {l, 0, 2 M}];
 xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, 2 M + 1}];
 Sp1[x_, q_] := 
  Sum[aS[i, j] pc[x, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aS1 pc1[x, q];
 Sp[x_] := 
  Sum[aS[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aS1 p1[x] + aS0;
 Ep1[x_, q_] := 
  Sum[aE[i, j] pc[x, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aE1 pc1[x, q];
 Ep[x_] := 
  Sum[aE[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aE1 p1[x] + aE0;
 Ip1[x_, q_] := 
  Sum[aI[i, j] pc[x, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aI1 pc1[x, q];
 Ip[x_] := 
  Sum[aI[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aI1 p1[x] + aI0;
 Ap1[x_, q_] := 
  Sum[aA[i, j] pc[x, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aA1 pc1[x, q];
 Ap[x_] := 
  Sum[aA[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aA1 p1[x] + aA0;
 Rp1[x_, q_] := 
  Sum[aR[i, j] pc[x, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aR1 pc1[x, q];
 Rp[x_] := 
  Sum[aR[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aR1 p1[x] + aR0;
 Mp1[x_, q_] := 
  Sum[aM[i, j] pc[x, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aM1 pc1[x, q];
 Mp[x_] := 
  Sum[aM[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aM1 p1[x] + aM0;
 varM = Join[aco0, aco1, 
   Flatten[Table[{aS[i, j], aE[i, j], aI[i, j], aA[i, j], aR[i, j], 
      aM[i, j]}, {j, 0, J, 1}, {i, 0, 2^j - 1, 1}]]];

 eq1[t_, q_] := -tn[q]/Gamma[1 - q] Sp1[t, q] + \[CapitalPi]p/
    Np0 - \[Mu]p Sp[t] - \[Eta]p Sp[
     t] (Ip[t] + \[Psi] Ap[t])/(Sp[t] + Ep[t] + Ip[t] + Ap[t] + 
       Rp[t]) - Np0 \[Eta]w Sp[t] Mp[t];
 eq2[t_, q_] := -tn[q]/Gamma[1 - q] Ep1[t, q] + \[Eta]p Sp[
     t] (Ip[t] + \[Psi] Ap[t])/(Sp[t] + Ep[t] + Ip[t] + Ap[t] + 
       Rp[t]) + 
   Np0 \[Eta]w Sp[t] Mp[t] - (1 - \[Theta]p) \[Omega]p Ep[
     t] - \[Theta]p \[Rho]p Ep[t] - \[Mu]p Ep[t];
 eq3[t_, q_] := -tn[q]/Gamma[1 - q] Ip1[t, 
     q] + (1 - \[Theta]p) \[Omega]p Ep[t] - (\[Tau]p + \[Mu]p) Ip[t];
 eq4[t_, q_] := -tn[q]/Gamma[1 - q] Ap1[t, q] + \[Theta]p \[Rho]p Ep[
     t] - (\[Tau]ap + \[Mu]p) Ap[t];
 eq5[t_, q_] := -tn[q]/Gamma[1 - q] Rp1[t, q] + \[Tau]p Ip[
     t] + \[Tau]ap Ap[t] - \[Mu]p Rp[t];
 eq6[t_, q_] := -tn[q]/Gamma[1 - q] Mp1[t, q] + \[CurlyRho]p Ip[
     t] + \[CurlyPi]p Ap[t] - \[Pi]p Mp[t];
 eq[q_] := 
  Flatten[Table[{eq1[t, q] == 0, eq2[t, q] == 0, eq3[t, q] == 0, 
     eq4[t, q] == 0, eq5[t, q] == 0, eq6[t, q] == 0}, {t, xcol}]];
 ParallelDo[ (* It is more efficient to parallelize this Do loop *)
  icv[i] = {Sp[0] == 8065518/Np0, Ep[0] == 200000/Np0, 
    Ip[0] == 282/Np0, Ap[0] == 200/Np0, Rp[0] == 0, 
    Mp[0] == 50000/Np0};
  eqM[i] = Join[eq[i], icv[i]];
  solv[i] = 
   FindRoot[eqM[i], Table[{varM[[j]], .1}, {j, Length[varM]}], 
    MaxIterations -> 1000];
  lstSv[i] = 
   Table[{x 120, Np0 Evaluate[Sp[x] /. solv[i]]}, {x, 0, 1, .01}];
  lstEv[i] = 
   Table[{x 120, Np0 Evaluate[Ep[x] /. solv[i]]}, {x, 0, 1, .01}];
  lstIv[i] = 
   Table[{x 120, Np0 Evaluate[Ip[x] /. solv[i]]}, {x, 0, 1, .01}];
  lstAv[i] = 
   Table[{x 120, Np0 Evaluate[Ap[x] /. solv[i]]}, {x, 0, 1, .01}];
  lstRv[i] = 
   Table[{x 120, Np0 Evaluate[Rp[x] /. solv[i]]}, {x, 0, 1, .01}];
  lstMv[i] = 
   Table[{x 120, Np0 Evaluate[Mp[x] /. solv[i]]}, {x, 0, 1, .01}];
  , {i, {99/100, 9/10, 8/10, 7/10, 6/10}}];
 (* Collect the definitions from the parallel kernels to the main \
one, e.g. by: (in principle this should work also by SharedFunction, \
but tat does not work well )*)
 Table[With[{lst = lst}, 
   DownValues[lst] = 
    Flatten@ParallelEvaluate[DownValues[lst]]], {lst, {lstSv, lstEv, 
    lstIv, lstAv, lstRv, lstMv}}];
 ]

Thank you! Actually I didn't see difference with my code if a piece of code in AbsoluteTiming[] using many times for parameters evaluation. In this case ParallelTable[] running even faster. — Alex Trounev, Commented May 18, 2020 at 10:27
It can be very tricky to optimize parallel evaluations in Mathematica. It also depends on the nature of the problem and also on hardware. — Rolf Mertig, Commented May 18, 2020 at 11:07
If memory serves, WaveletPsi[] and WaveletPhi[] are supposed to produce Piecewise[] functions, so if need be, someone could convert the function (by hand?) into something with Which[] (or re-express in terms of UnitStep[]) so that it could be compileable. — J. M.'s missing motivation, Commented May 18, 2020 at 11:58
@Alex, that must have been a recent thing (I'm still on version 11); if it is working for you, then my concerns are moot. But just to check: you can use CompilePrint[] on your functions and check for MainEvaluate calls in the output. — J. M.'s missing motivation, Commented May 18, 2020 at 17:01

Stack Exchange Network

Solver for COVID-19 epidemic model with the Caputo fractional derivatives

1 Answer 1

Your Answer

Not the answer you're looking for? Browse other questions tagged
differential-equations
calculus-and-analysis
or ask your own question.

Linked

Hot Network Questions

Solver for COVID-19 epidemic model with the Caputo fractional derivatives

1 Answer 1

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged differential-equationscalculus-and-analysis or ask your own question.

Linked

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
differential-equations
calculus-and-analysis
or ask your own question.