Nested For loop to plot multiple plots on a single graph

Please help me to plot these 8 graphs on the same figure by changing the values of eta and epsilon. I can plot it separately and use Show to draw it together, but I wanted to use loop to get the all cases on a single figure using hold on and hold off but it seems not to work.

    s = 2.6 10^4;
    d = 0.0026;
   \[Beta] = 3 10^(-7);
    \[Delta] = 0.5;
    c = 5;
    p = 100;
    \[Eta] = {0.00, 1.00};
    \[Epsilon] = {0.80, 0.95, 0.99, 1};
    eqn1 = {0, 0, 0, 0};
    (*Initial Conditions*)
    t0 = (c*\[Delta])/(p*\[Beta]); x0 = 2.0*  10^(6); v0 = 10^7;
    For[i = 1, j = 1; i < 3, j < 5; j++, i++; 
     eqn1[[i]] := {T'[t] == 
        s - d * T[t] - (1 - \[Eta][[i]])*\[Beta] *T [t]*V[t], 
       x'[t] == (1 - \[Eta][[i]])*\[Beta] *T[t]*V[t] - \[Delta] *x[t], 
       V'[t] == (1 - \[Epsilon][[j]])*p *x[t] - c *V[t], T[0] == t0, 
       x[0] == x0, V[0] == v0}; 
     soln11[[i]] = NDSolve[eqn1[[i]], {T, x, V}, {t, 0, 14}]; 
     p1 = LogPlot[
       Evaluate[{V[t]} /. First@soln11[[i]], {t, 0, 14}, 
        PlotRange -> {100, 10^8}, PlotPoints -> 200, Frame -> True, 
        PlotStyle -> {Orange, Thick}]]

neumannwidloop.nb (7.33 KB)

s = 2.6 10^4;(cells per ml per day))
d = 0.0026;(day^(-1))
[Beta] = 3 10^(-7);((Virion per mililiter)^(-1))
[Delta] = 0.5;(day^(-1))
c = 5;(day^(-1))
p = 100;(virions per ml per cell per day)
[Eta] = {0.00, 1.00};
[Epsilon] = {0.80, 0.95, 0.99, 1};
eqn1 = {0, 0, 0, 0};
(Initial Conditions)
t0 = (c*[Delta])/(p*[Beta]); x0 = 2.0* 10^(6); v0 = 10^7;
For[i = 1, j = 1; i < 3, j < 5; j++, i++;
eqn1[[i]] := {T’[t] ==
s - d * T[t] - (1 - [Eta][[i]])*[Beta] *T [t]V[t],
x’[t] == (1 - [Eta][[i]])
[Beta] *T[t]*V[t] - [Delta] *x[t],
V’[t] == (1 - [Epsilon][[j]])*p *x[t] - c *V[t], T[0] == t0,
x[0] == x0, V[0] == v0};
soln11[[i]] = NDSolve[eqn1[[i]], {T, x, V}, {t, 0, 14}];
p1 = LogPlot[
Evaluate[{V[t]} /. First@soln11[[i]], {t, 0, 14},
PlotRange → {100, 10^8}, PlotPoints → 200, Frame → True,
PlotStyle → {Orange, Thick}]]

For reasons unknown to me the above code didn’t work.
I do not understand why the i-loop starts with 2 and j doesn’t reach 4:

For[i = 1, j = 1; i < 3, j < 5; j++, i++;
  Print[" i ", i];
  Print[" j ", j]
  ];

Here it does

For[i = 1, i < 3, i++,
 For[j = 1, j < 5, j++,
  Print[" i ", i];
  Print[" j ", j]
  ]
 ]

But what about this (without For…)?

s = 2.6 10^4;(*cells per ml per day)*)
d = 0.0026;(*day^(-1)*)
\[Beta] =  3 10^(-7);(*(Virion per mililiter)^(-1)*)
\[Delta] =0.5;(*day^(-1)*)
c = 5;(*day^(-1)*)
p = 100;(*virions per ml per cell per day*)
\[Eta] = {0.00, 1.00};
\[Epsilon] = {0.80, 0.95, 0.99, 1};
       (*Initial Conditions*)
t0 = (c*\[Delta])/(p*\[Beta]); x0 = 2.0*10^(6); v0 = 10^7;

eqn1a[i_,  j_] := {T'[t] == s - d*T[t] - (1 - \[Eta][[i]])*\[Beta]*T[t]*V[t], 
  x'[t] == (1 - \[Eta][[i]])*\[Beta]*T[t]*V[t] - \[Delta]*x[t], 
  V'[t] == (1 - \[Epsilon][[j]])*p*x[t] - c*V[t], T[0] == t0, 
  x[0] == x0, V[0] == v0}

soln11[eq_] := NDSolve[eq, {T, x, V}, {t, 0, 14}]

tt = Flatten[Table[eqn1a[i, j], {i, 2}, {j, 4}], 1];
sol = Flatten[soln11 /@ tt];
solv = (V /. #)[t] & /@ Select[sol, MemberQ[#, V] &];

LogPlot[solv, {t, 0, 14}, PlotRange -> {100, 10^8}, PlotPoints -> 200,
  Frame -> True, PlotStyle -> {Orange, Thick}]

Similarly

solx = (x /. #)[t] & /@ Select[sol, MemberQ[#, x] &];    
solT = (T /. #)[t] & /@ Select[sol, MemberQ[#, T] &];

Thank you Sir!! It seems working.

Hi Hans

I do not understand why the i-loop starts with 2 and j doesn’t reach 4.

In WL the role of ; and , is reversed relative to languages like C.

For[i = 1; j = 1, i < 3, j++; i++,
  Print[" i ", i];
  Print[" j ", j]];

As you suggested, Table is a much better option that For.

“In WL the role of ; and , is reversed relative to languages like C.”

OK, but

For[i = 1; j = 1, i < 3, j++; i++, Print[" i ", i];
  Print[" j ", j]];


 i 1    
 j 1    
 i 2    
 j 2

Obviously you miss the pairs i = 1, j = 2 and so on, and j never reaches 4. But I think Abiy wanted all combinations of i and j. Seems this construct isn’t appropriate for what he intended.

Look at

For[i = 1; j = 1, i < 3; j < 5, j++; i++, Print[{i, j}]];

{1,1}

{2,2}

{3,3}

{4,4}

which is obviously equivalent to

For[i = 1; j = 1, Or[i < 3, j < 5], j++; i++, Print[{i, j}]];

{1,1}

{2,2}

{3,3}

{4,4}

then, and this is according to what is said in the Help-Section : The sequence of evaluation is test, body, incr .

So
i and j are set to 1, the test gives True , and i AND j are incremented and so on, until the test yields False. Try setting j to 3 and substitute j<5 with j < 10 for example.

I want to have a dashed lines for eta=1 and thick for eta=0 like on the code below. But without the for loop, it becomes redundant. So Can you please suggest a better way?

  s = 2.6 10^4;
d = 0.0026;
\[Beta] =3 *10^(-7);\[Delta] = \
        0.5;c = 5;p = 100;\[Eta] = {0.00, 1.00};
        \[Epsilon] = {0.80, 0.95, 0.99, 1};
        (*Initial Conditions*)
        t0 = (c*\[Delta])/(p*\[Beta]); x0 = 2.0*10^(6); v0 = 10^7;
        
        eqn1a[i_, 
          j_] := {T'[t] == s - d*T[t] - (1 - \[Eta][[i]])*\[Beta]*T[t]*V[t], 
          x'[t] == (1 - \[Eta][[i]])*\[Beta]*T[t]*V[t] - \[Delta]*x[t], 
          V'[t] == (1 - \[Epsilon][[j]])*p*x[t] - c*V[t], T[0] == t0, 
          x[0] == x0, V[0] == v0}
        
        soln11[eq_] := NDSolve[eq, {T, x, V}, {t, 0, 14}]
        
        tt1 = Flatten[Table[eqn1a[i, j], {i, 1}, {j, 4}], 1];
        sol1 = Flatten[soln11 /@ tt1];
        solv1 = (V /. #)[t] & /@ Select[sol1, MemberQ[#, V] &];
        
        pfig1 = LogPlot[solv1, {t, 0, 14}, PlotRange -> {3*100, 2*10^7}, 
           PlotPoints -> 200, Frame -> True, PlotStyle -> {Orange, Thick}];
        soln12[eq_] := NDSolve[eq, {T, x, V}, {t, 0, 14}]
        
        tt2 = Flatten[Table[eqn1a[i, j], {i = 2}, {j, 4}], 1];
        sol2 = Flatten[soln12 /@ tt2];
        solv2 = (V /. #)[t] & /@ Select[sol2, MemberQ[#, V] &];
        
        pfig2 = LogPlot[solv2, {t, 0, 14}, PlotRange -> {3*100, 2*10^7}, 
           PlotPoints -> 200, Frame -> True, PlotStyle -> {Blue, Dashed}];
        Show[pfig1, pfig2]

Do you mean this? (Note that the first four functions are for eta = 0 )

s = 2.6 10^4;(*cells per ml per day)*)d = 0.0026;(*day^(-1)*)\[Beta] =
  3 10^(-7);(*(Virion per mililiter)^(-1)*)\[Delta] = \
0.5;(*day^(-1)*)c = 5;(*day^(-1)*)p = 100;(*virions per ml per cell \
per day*)\[Eta] = {0.00, 1.00};
\[Epsilon] = {0.80, 0.95, 0.99, 1};
(*Initial Conditions*)
t0 = (c*\[Delta])/(p*\[Beta]); x0 = 2.0*10^(6); v0 = 10^7;

eqn1a[i_, 
  j_] := {T'[t] == s - d*T[t] - (1 - \[Eta][[i]])*\[Beta]*T[t]*V[t], 
  x'[t] == (1 - \[Eta][[i]])*\[Beta]*T[t]*V[t] - \[Delta]*x[t], 
  V'[t] == (1 - \[Epsilon][[j]])*p*x[t] - c*V[t], T[0] == t0, 
  x[0] == x0, V[0] == v0}

soln11[eq_] := NDSolve[eq, {T, x, V}, {t, 0, 14}]

tt = Flatten[Table[eqn1a[i, j], {i, 2}, {j, 4}], 1];
sol = Flatten[soln11 /@ tt];
solv = (V /. #)[t] & /@ Select[sol, MemberQ[#, V] &];

p1 = LogPlot[solv, {t, 0, 14}, PlotRange -> {100, 10^7}, 
  PlotPoints -> 200, Frame -> True, 
  PlotStyle -> Table[If[i <= 4, Thick, Dashed], {i, 8}]]

Yes Sir you understood me well, but the above code makes all 8 plots dashed, no thick line (just changed the thick lines of the first code by dashed, I want it to be mixed: first four thick and the rest Dashed).

In fact, I am afraid that I don’t understand you. What do you mean?
And what do you say considering (in my code)

p1 = LogPlot[solv, {t, 0, 14}, PlotRange -> {100, 10^7}, 
  PlotPoints -> 200, Frame -> True, 
  PlotStyle -> Table[If[ i <= 4, {Red, Thickness[.01]}, {Blue, Thickness[.008], Dashed}], {i, 8}]]

You might as well want to look at

first four functions:

p2 = LogPlot[Evaluate[Join[Take[solv, 4], {1, 1, 1, 1}]], {t, 0, 14},
PlotRange -> {100, 10^7}, PlotPoints -> 200, Frame -> True, 
PlotStyle -> Table[If[ i <= 4, {Red, Thickness[.01]}, {Blue, Thickness[.005],  Dashed}], {i, 8}]]

last four functions

p3 = LogPlot[Evaluate[Join[{1, 1, 1, 1}, Take[solv, -4]]], {t, 0, 14},
   PlotRange -> {100, 10^7}, PlotPoints -> 200, Frame -> True, 
  PlotStyle -> Table[If[ i <= 4, {Red, Thickness[.01]}, {Blue, Thickness[.008],  Dashed}], {i, 8}]]

I would like to plot V and VNI together as V+VNI and the code below is not helping. Mr. Hans Dolhaine, May I get the usual help please?

  s = 2.6 10^4;
      d = 0.0026;
      \[Beta] = 2.25 10^(-7);
      \[Delta] = 1;
      c = 6;
      p = 2.9;  
    \[Epsilon] = {0.96, 0.5};
        \[Rho] = {0, 1};
        eqn1g[i_, j_] := {x'[t] == \[Beta] *T*V[t] - \[Delta] *x[t], 
           V'[t] == (1 - \[Rho][[i]]) (1 - \[Epsilon][[j]])*p *x[t] - 
             c *V[t], 
           VNI'[t] == \[Rho][[i]]*(1 - \[Epsilon][[j]])*p *x[t] - c *VNI[t], 
           x[0] == x0, V[0] == v0, VNI[0] == 0};
        
        soln112[eq_] := NDSolve[eq, {x, V, VNI}, {t, 0, 14}];
        tableway1 = Flatten[Table[eqn1g[i, j], {i, 2}, {j, 2}], 1];
        Solsn1 = Flatten[soln112 /@ tableway];
        Solv1 = (V /. #)[t] & /@ Select[Solsn1, MemberQ[#, V] &];
        Solv2 = (VNI /. #)[t] & /@ Select[Solsn1, MemberQ[#, VNI] &];
        p3 = LogPlot[
          Evaluate[Join[Take[Solv1 + Solv2, 4], {1, 1, 1, 1}]], {t, 0, 14}, 
          PlotRange -> {100, 10^7}, PlotPoints -> 200, Frame -> True, 
          FrameLabel -> {"Time" , "copies "}]

Good evening Abiy,

You have a typo in your code:

Solsn1 = Flatten[soln112 /@ tableway];

there you should refer to

tableway1

as defined further above. Moreover you have to define values for v0, x0 and T. I guess T means Temperature and I assigned a value of 273.

With these modifications I arrive at a code which gives a plot. Is it this what you meant?

s = 2.6 10^4;
d = 0.0026;
\[Beta] = 2.25 10^(-7);
\[Delta] = 1;
c = 6;
p = 2.9;
\[Epsilon] = {0.96, 0.5};
\[Rho] = {0, 1};
v0 = 10^7;
x0 = 2. 10^6;
T = 273;
eqn1g[i_, j_] := {x'[t] == \[Beta]*T*V[t] - \[Delta]*x[t], 
   V'[t] == (1 - \[Rho][[i]]) (1 - \[Epsilon][[j]])*p*x[t] - c*V[t], 
   VNI'[t] == \[Rho][[i]]*(1 - \[Epsilon][[j]])*p*x[t] - c*VNI[t], 
   x[0] == x0, V[0] == v0, VNI[0] == 0};

soln112[eq_] := NDSolve[eq, {x, V, VNI}, {t, 0, 14}];

tableway1 = Flatten[Table[eqn1g[i, j], {i, 2}, {j, 2}], 1];
Solsn1 = Flatten[soln112 /@ tableway1];

Solv1 = (V /. #)[t] & /@ Select[Solsn1, MemberQ[#, V] &];
Solv2 = (VNI /. #)[t] & /@ Select[Solsn1, MemberQ[#, VNI] &];
p3 = LogPlot[
  Evaluate[Join[Take[Solv1 + Solv2, 4], {1, 1, 1, 1}]], {t, 0, 14}, 
  PlotRange -> {100, 10^7}, PlotPoints -> 200, Frame -> True, 
  FrameLabel -> {"Time", "copies "}]

Sure correcting the typo error solved everything. Thank you for you quick and unreserved help. U are helping me get better with Mathematica everyday. Big respect SIR!!

I want to have the y-axis ticks at {0,0.2,0.4,0.6,0.8,1, 1.2,1.4, 1.6} instead of in {0,10^(-4), 0.001, 0.01, 0.1,1}. Please see the plot of the code below and help me in improving it by plotting the y-axis in step size of 0.2. Both Ticks and Frame Ticks were not able to help me.

 k = 6.7 *10^-8; c = 
         2*10^6; \[Alpha] = 1.34; \[Lambda] = 3.3; \[Mu] = 2; X = 10^6;
        r = 0.17*10^6;
        S = 0.83 *10^6;
        P = 10^7;
        a = \[Alpha]/(k*c);
         l = \[Lambda]/(k*c);
        m = \[Mu]/(k*c);
        s1 = S/c; R1 = r/c; p1 = P/c;
        b = 15;
    (*Initial Conditions*)
    
    x0 = 0.3; y0 = 0.2; z0 = 5;
    
    (*Differential Equations*)
    
    rhs1[x_, y_, z_] := x (a (1 - (x + y)) - z);
    rhs2[x_, y_, z_] := -l y + x z;
    rhs3[x_, y_, z_] := -z x + b l y - m z;

   detstab := {x'[t] == rhs1[x[t], y[t], z[t]], 
       y'[t] == rhs2[x[t], y[t], z[t]], z'[t] == rhs3[x[t], y[t], z[t]], 
       x[0] == x0, y[0] == y0, z[0] == z0};
    soln11 = NDSolve[detstab, {x, y, z}, {t, 0, 10}];
    LogPlot[Evaluate[{x[t], y[t], z[t]} /. First@soln11], {t, 0, 10}, 
     PlotRange -> {{0, 10}, {0, 1.2}}, 
     Ticks -> {Table[i + 2, {i, 0, 10}], Table[i + 0.2, {i, 0, 1}]}, 
     Frame -> True]

Hi Abiy,

In a LogPlot the y axis is always logarithmic. You could plot the Log of the function on a linear y scale. e.g.

Plot[Evaluate[Log@{x[t], y[t], z[t]} /. First@soln11], {t, 0, 10}, 
 Frame -> True, PlotLegends -> {"x", "y", "z"}]

Or a linear plot with a specific range / ticks

Plot[Evaluate[{x[t], y[t], z[t]} /. First@soln11], {t, 0, 10},
 PlotRange -> {0, 1.6},
 Frame -> True,
 FrameTicks -> {{Range[0, 1.6, .2], None}, {Automatic, None}},
 PlotLegends -> {"x", "y", "z"}]

Dear Rohit Namjoshi, thank you, this is what I wanted. However, I have one more question, as you can see below I have been using LogPlot and the y-axis keeps giving me logarithmic values but what I need is to plot the Log of the function on a linear y scale, as your first suggestion, I tried to improve the code but keeps giving me a wrong result. So will you help in changing the codes below to plot the Log of the function instead of using Logplot.

s = 2.6 10^4;
d = 0.0026;
\[Beta] = 2.25 10^(-7);
\[Delta] = 1;
c = 6;
p = 2.9;
\[Epsilon] = {0.96, 0.5};
\[Rho] = {0, 1};
v0 = 10^7;
x0 = 2. 10^6;
T = 273;
eqn1g[i_, j_] := {x'[t] == \[Beta]*T*V[t] - \[Delta]*x[t], 
   V'[t] == (1 - \[Rho][[i]]) (1 - \[Epsilon][[j]])*p*x[t] - c*V[t], 
   VNI'[t] == \[Rho][[i]]*(1 - \[Epsilon][[j]])*p*x[t] - c*VNI[t], 
   x[0] == x0, V[0] == v0, VNI[0] == 0};

soln112[eq_] := NDSolve[eq, {x, V, VNI}, {t, 0, 14}];

tableway1 = Flatten[Table[eqn1g[i, j], {i, 2}, {j, 2}], 1];
Solsn1 = Flatten[soln112 /@ tableway1];

Solv1 = (V /. #)[t] & /@ Select[Solsn1, MemberQ[#, V] &];
Solv2 = (VNI /. #)[t] & /@ Select[Solsn1, MemberQ[#, VNI] &];
p3 = LogPlot[
  Evaluate[Join[Take[Solv1 + Solv2, 4], {1, 1, 1, 1}]], {t, 0, 14}, 
  PlotRange -> {100, 10^7}, PlotPoints -> 200, Frame -> True, 
  FrameLabel -> {"Time", "copies "}]

Hi Abiy,

What is the reason for joining {1, 1, 1, 1} to the list of plots? It is just going to plot a straight line at y = 1. Is this what you are looking for?

p3 = Plot[Evaluate[Log@Take[Solv1 + Solv2, 4]], {t, 0, 14},
  PlotRange -> Log@{100, 10^7},
  PlotPoints -> 200,
  Frame -> True,
  FrameLabel -> {"Time", "copies "}]

I appreciate the help. Sure joining is not important !!

@Rohit: I introduced

Join[Take[solv, 4], {1, 1, 1, 1}]

in a post far above to have easy acces to the “dashed lines” by then using (see p3 far above)

Join[{1,1,1,1},Take[solv, - 4]]

avoiding to have to change (which of course could be done as well)

PlotStyle -> Table[If[ i <= 4, {Red, Thickness[.01]}, {Blue, Thickness[.008], Dashed}], {i, 8}]]

and knowing that Log[1] == 0 will coincide with the x-axis and not disturb the plot.

Abyi wanted to have thick and dashed lines.

It seems that this is not necessary later on.

@Abiy : By the way, you said

I would like to plot V and VNI together as V+VNI

Look at your equations and initial values (with the input in the code above):

In[40]:= eee = Flatten[Table[eqn1g[i, j], {i, 2}, {j, 2}]];
Select[eee, MemberQ[#, VNI'[t]] &]
Select[eee, MemberQ[#, VNI[0]] &]

Out[41]= {
\!\(\*SuperscriptBox["VNI", "\[Prime]",
MultilineFunction->None]\)[t] == -6 VNI[t], 
\!\(\*SuperscriptBox["VNI", "\[Prime]",
MultilineFunction->None]\)[t] == -6 VNI[t], 
\!\(\*SuperscriptBox["VNI", "\[Prime]",
MultilineFunction->None]\)[t] == -6 VNI[t] + 0.116 x[t], 
\!\(\*SuperscriptBox["VNI", "\[Prime]",
MultilineFunction->None]\)[t] == -6 VNI[t] + 1.45 x[t]}

Out[42]= {VNI[0] == 0, VNI[0] == 0, VNI[0] == 0, VNI[0] == 0}

The first two VNI can be integrated immediately giving

VNI[t_] := VNI[0] Exp[- 6 t]

and with

VNI[0] == 0

this is always zero (check it with the results of NDSolve).

A similar direct integration is possible for the last two V’s

Dear Hans Dolhaine, you are right that I needed it at the beginning not to repeat the same code for different plots( dashed(p2) and thick(p3)) using an if loop(which can easily be applied on Matlab). But the if loop were not able to give me the plot I want(had the if loop were working p2 would have been enough to plot all and p3 was not needed). Thus I finally decided to ignore the join as Rohit suggested and plot thick and dashed lines on different codes as below. I honestly thank both of you for being a big help in my Mathematica coding. If possible, help me to keep p2 and p3 in a single code (i=1(j=1:4)thick lines and i=2(j=1:4) dashed lines) using if loop than being redundant and not to use the Show function. Here is the code for both to see and if possible improve.
Kind regards. (Thank you again!!)

s = 2.6*10^4;
d = 0.0026;
\[Beta] = 3 *10^(-7);
\[Delta] = {0.5, 5};
c = 5;
p = {5, 100};

\[Epsilon] = {0.00, 0.00};
\[Eta] = {0.80, 0.95, 0.99, 1};

\[Eta]a = {0.00, 1.00};
\[Epsilon]a = {0.80, 0.95, 0.99, 1};
(*Initial Conditions*)
t0 = {(c*\[Delta][[2]])/(p[[2]]*\[Beta]), (c*\[Delta][[1]])/(p[[
       2]]*\[Beta])};
x0 = 2.0*10^(6); v0 = 4.0*10^7;
(* SDES for p=100, delta=0.5 for fig A*)
eqn1a[i_, 
  j_] := {T'[t] == s - d*T[t] - (1 - \[Eta]a[[i]])*\[Beta]*T[t]*V[t], 
  x'[t] == (1 - \[Eta]a[[i]])*\[Beta]*T[t]*V[t] - \[Delta][[1]]*x[t], 
  V'[t] == (1 - \[Epsilon]a[[j]])*p[[2]]*x[t] - c*V[t], 
  T[0] == t0[[2]], x[0] == x0, V[0] == v0}

 soln11[eq_] := NDSolve[eq, {T, x, V}, {t, 0, 14}]

(* Solving the SDES *) 
tt2 = Flatten[Table[eqn1a[i, j], {i, 2}, {j, 4}], 1];
sol2 = Flatten[soln11 /@ tt2];
solv2 = (V /. #)[t] & /@ Select[sol2, MemberQ[#, V] &];

p2 = Plot[Evaluate[Log10@Take[solv2, 4]], {t, 0, 14}, 
   PlotRange -> {2.6, 8}, Frame -> True, 
   FrameLabel -> {Days , "Simulated Decrease Log10 HCV"},
   PlotStyle -> Table[If[i <= 4, {Red, Thin}, {Blue, "--"}], {i, 8}]];
(*Plot B eta=0 *)
p3 = Plot[Evaluate[Log10@Take[solv2, -4]], {t, 0, 14}, 
   PlotRange -> {2.6, 8}, Frame -> True, 
   FrameLabel -> {Days , "Simulated Decrease Log10 HCV"}, 
   PlotStyle -> Table[If[i <= 4, {Blue, Dashed}], {i, 8}]];
(* Draw together*)
p41 = Show[p2, p3]