Try

Table[Which[i == 1 \[Or] i == 3 \[Or] i == 5, {Red, Thin}, 
  i == 2 \[Or] i == 4 \[Or] i == 6, {Green, Thin}, 
  True, {Blue, Thin}], {i, 8}]

Dear Hans Dolhaine, my version is 11 and it has a way of solving stochastic differential equations!! ItoProcess is also known as Ito diffusion or stochastic differential equation (SDE). ItoProcess is a continuous-time and continuous-state random process. If the drift a is an -dimensional vector and the diffusion b an × -dimensional matrix, the process is -dimensional and driven by an -dimensional WienerProcess.

An Ito process is a type of stochastic process described by Japanese mathematician Kiyoshi Itô, which can be written as the sum of the integral of a process over time and of another process over a Brownian motion.

My problem is I keep getting this error, "Thread::tdlen: Objects of unequal length in 1 cannot be combined." This message is generated when a threading operation is used with expressions that do not have the same length. I just am trying to find a way to run it.

But thank you for your help till now.

Dear Hans Dolhaine and Rahul, I need one more favor as always!! I wanted to change the equation (eqn1a) as ItoProcess so that I can change the above equations by adding some stochastic terms. Then instead of NDSolve, I want to use Randomfunction and finally draw the ListLinePlot of solv2. But it seems something is wrong and I cannot get any results. Can you please run the code below and lead me in the right direction??? Thank you!!

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*)
t01 = 2.4*10^6;
v01 = 4.0*10^5;
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*)
Subscript[\[Sigma], 1] = 0.07; Subscript[\[Sigma], 2] = 0.5; \
Subscript[\[Sigma], 3] = 3;


s2 = s/d; R2 = 0; p2 = 0;

eqn1a[i_, j_] := 
 ItoProcess[{\[DifferentialD]T[
      t] == (s - 
        d*T[t] - (1 - \[Eta]a[[i]])*\[Beta]*T[t]*
         V[t]) \[DifferentialD]t + 
     Subscript[\[Sigma], 
      1] (x[t] - s2) \[DifferentialD]w1[t], \[DifferentialD]x[
      t] == ((1 - \[Eta]a[[i]])*\[Beta]*T[t]*V[t] - \[Delta][[1]]*
         x[t]) \[DifferentialD]t + 
     Subscript[\[Sigma], 
      2] (y[t] - R2) \[DifferentialD]w2[t], \[DifferentialD]V[
      t] == ((1 - \[Epsilon]a[[j]])*p[[2]]*x[t] - 
        c*V[t]) \[DifferentialD]t + 
     Subscript[\[Sigma], 3] (z[t] - p2) \[DifferentialD]w3[t]}, {x[t],
    y[t], z[t]}, {{x, y, z}, {t0[[2]], x0, v0}}, 
  t, {w1 \[Distributed] WienerProcess[], 
   w2 \[Distributed] WienerProcess[] , 
   w3 \[Distributed] WienerProcess[]}]

soln11[eq_] := 
  RandomFunction[eq, {T, x, V}, {t, 0., 14., 0.0001}, 
   Method -> "Milstein"];
tt2 = Flatten[Table[eqn1a[i, j], {i, 2}, {j, 4}], 1];
sol2 = Flatten[soln11 /@ tt2];
solv2 = (V /. #)[t] & /@ Select[sol2, MemberQ[#, V] &];


p2 = ListLinePlot[Evaluate[solv2], PlotRange -> All, Frame -> True, 
  FrameLabel -> {Days, "Simulated Decrease Log10 HCV"}, 
  PlotStyle -> Table[If[i <= 4, {Red, Thin}, {Blue, Dashed}], {i, 8}]]

Dear Abiy,

if I understand you correctly you want a single plot-command.

Try a replacement of your p2 in the above code by (note that I dropped the ; at the end)

p2 = Plot[Evaluate[Log10 /@ solv2], {t, 0, 14}, PlotRange -> {2.6, 8},
   Frame -> True, 
  FrameLabel -> {Days, "Simulated Decrease Log10 HCV"}, 
  PlotStyle -> Table[If[i <= 4, {Red, Thin}, {Blue, Dashed}], {i, 8}]]

and delete p3 and Show[...]

(We had this already far above with p1 but as LogPlot and differring thicknesses and dashing-style)

solv2 contains all the information for your plot (8 functions). You don't need to pick the first four, plot them, then pick the last four, plot them and unite the plots with Show.

@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

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

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

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

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).

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]

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.

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] &];