(* Wolfram Mathematica Code *)

CoinFlip[p_, n_: 1] := Sum[If[RandomReal[] < p, 1, 0], {n}]
CreateHousehold[hp_: 2, hc_: 3, wealth_: 500, sick_: False] :=
 {1 + CoinFlip[(hp - 1)/3, 3], 0, 0, 0, CoinFlip[hc/7, 7], 0, 0, 0, 
   wealth} + If[sick, {-1, 1, 0, 0, 0, 0, 0, 0, 0}, 0]
Monte[Tp_, Tc_, CC_: 30, CT_: 50, PP_: 50, wealth_: 500, tmax_: 120] :=
  Module[
   {
    n = 500, hp = 2, hc = 3, (* Formerly arguments *)
    \[Beta]p = 0.3/2500,
    \[Beta]c = 0.5/2500,
    \[Gamma]p = 0.1,
    \[Gamma]c = 0.1,
    \[Gamma]pt = 0.5,
    \[Gamma]ct = 0.5,
    \[Iota]p = 0.03/2500,
    \[Iota]c = 0.05/2500,
    \[Beta]p2 = 0.4/2500,
    \[Beta]c2 = 0.7/2500,
    \[Iota]p2 = 0.04/2500,
    \[Iota]c2 = 0.07/2500,
    Sx = {1, 5},
    Ix = {2, 3, 6, 7},
    Rx = {4, 8},
    Px = {1, 2, 3, 4},
    Cx = {5, 6, 7, 8},
    Mx = {9},
    Population, t, newpop, Iall, hh, ip, ic, ipt, ict, rp, rpt, rc, 
    rct, i,
    MONEY
    },
   Population[0] = 
    Append[Table[CreateHousehold[hp, hc, wealth], {n - 1}], 
     CreateHousehold[hp, hc, wealth, True]];
   For[t = 1, t <= tmax, t++,
    Iall[t - 1] = Plus @@ Plus @@ Population[t - 1][[All, Ix]];
    newpop = Population[t - 1];
    For[i = 1, i <= n, i++,
     hh = Population[t - 1][[i]];
     ip = 
      CoinFlip[If[hh[[9]] > 0, \[Beta]p, \[Beta]p2] Iall[t - 1], 
       hh[[1]]];
     ic = 
      CoinFlip[If[hh[[9]] > 0, \[Beta]c, \[Beta]c2] Iall[t - 1], 
       hh[[5]]];
     ip += 
      CoinFlip[If[hh[[9]] > 0, \[Iota]p, \[Iota]p2] Plus @@ hh[[Ix]], 
       hh[[1]]];
     ic += 
      CoinFlip[If[hh[[9]] > 0, \[Iota]c, \[Iota]c2] Plus @@ hh[[Ix]], 
       hh[[5]]];
     ip = Min[ip, hh[[1]]];
     ic = Min[ic, hh[[1]]];
     ipt = CoinFlip[Tp, ip];
     ict = CoinFlip[Tc, ic];
      If[hh[[9]] < 0 && CT > 0, (* This cut-off was changed to zero, 
      instead of CT(ipt+ict) *)
      ipt = 0; ict = 0;
      ];
     ic -= ict;
     ip -= ipt;
     newpop[[i]] += {-ip - ipt, ip, ipt, 0, -ic - ict, ic, ict, 0, 0};
     rp = CoinFlip[\[Gamma]p, hh[[2]]];
     rpt = CoinFlip[\[Gamma]pt, hh[[3]]];
     rc = CoinFlip[\[Gamma]c, hh[[6]]];
     rct = CoinFlip[\[Gamma]ct, hh[[7]]];
     newpop[[i]] += {0, -rp, -rpt, rp + rpt, 0, -rc, -rct, rc + rct, 
       0};
     newpop[[i, 9]] += 
      hh.{PP, PP - CT, -CC, PP, -CC, -CC - CT, -CC, -CC, 0};
     ];
    Population[t] = newpop;
    ];
   Return[
    Table[Population[t], {t, 0, tmax}]
    ];
   ];
MyODE[Tp_, Tc_, CC_: 30, CT_: 50, PP_: 50, wealth_: 500, tmax_: 120] :=
  Module[{
    n = 500, hp = 2, hc = 3,(* formerly arguments *)
    ODEDATA, h,
    \[Beta]p, \[Beta]c, \[Gamma]p, \[Gamma]c, \[Gamma]pt, \[Gamma]ct, \
\[Iota]p, \[Iota]c, \[Beta]p2, \[Beta]c2, \[Iota]p2, \[Iota]c2,
    Sp, Sc, Ip, Ic, Ipt, Ict, Rp, Rc, Mh, IpC, IcC, IptC, IctC
    },
   h = hc + hp;
   \[Gamma]p = -Log[0.9];
   \[Gamma]c = -Log[0.9];
   \[Gamma]pt = -Log[0.5];
   \[Gamma]ct = -Log[0.5];
   \[Beta]p = \[Gamma]p/0.1 0.3/(n h);
   \[Beta]c = \[Gamma]p/0.1 0.5/(n h);
   \[Iota]p = 0.03/(n h);
   \[Iota]c = 0.05/(n h);
   \[Beta]p2 = 0.4/(n h);
   \[Beta]c2 = 0.7/(n h);
   \[Iota]p2 = 0.04/(n h);
   \[Iota]c2 = 0.07/(n h);
   ODEDATA = NDSolveValue[
     {
      Sp'[
        t] == -Piecewise[{{\[Beta]p, Mh[t] > 0}, {\[Beta]p2, 
             Mh[t] <= 0}}] Sp[t] (Ip[t] + Ipt[t] + Ic[t] + Ict[t]) - 
        Piecewise[{{\[Iota]p, Mh[t] > 0}, {\[Iota]p2, Mh[t] <= 0}}] ( 
          h/(n h - h)) Sp[t] (Ip[t] + Ipt[t] + Ic[t] + Ict[t]),
      Sc'[
        t] == -Piecewise[{{\[Beta]c, Mh[t] > 0}, {\[Beta]c2, 
             Mh[t] <= 0}}] Sc[t] (Ip[t] + Ipt[t] + Ic[t] + Ict[t]) - 
        Piecewise[{{\[Iota]c, Mh[t] > 0}, {\[Iota]c2, Mh[t] <= 0}}] ( 
          h/(n h - h)) Sc[t] (Ip[t] + Ipt[t] + Ic[t] + Ict[t]),
      Ip'[t] == 
       Piecewise[{{1 - Tp, CT == 0 || Mh[t] > 0}, {1, 
            Mh[t] <= 
             0}}] (Piecewise[{{\[Beta]p, Mh[t] > 0}, {\[Beta]p2, 
               Mh[t] <= 0}}] Sp[t] (Ip[t] + Ipt[t] + Ic[t] + Ict[t]) +
            Piecewise[{{\[Iota]p, Mh[t] > 0}, {\[Iota]p2, 
               Mh[t] <= 0}}] ( h/(n h - h)) Sp[
             t] (Ip[t] + Ipt[t] + Ic[t] + Ict[t])) - \[Gamma]p Ip[t],
      Ic'[t] == 
       Piecewise[{{1 - Tc, CT == 0 || Mh[t] > 0}, {1, 
            Mh[t] <= 
             0}}] (Piecewise[{{\[Beta]c, Mh[t] > 0}, {\[Beta]c2, 
               Mh[t] <= 0}}] Sc[t] (Ip[t] + Ipt[t] + Ic[t] + Ict[t]) +
            Piecewise[{{\[Iota]c, Mh[t] > 0}, {\[Iota]c2, 
               Mh[t] <= 0}}] ( h/(n h - h)) Sc[
             t] (Ip[t] + Ipt[t] + Ic[t] + Ict[t])) - \[Gamma]c Ic[t],
      Ipt'[t] == 
       Piecewise[{{Tp, CT == 0 || Mh[t] > 0}, {0, 
            Mh[t] <= 
             0}}] (Piecewise[{{\[Beta]p, Mh[t] > 0}, {\[Beta]p2, 
               Mh[t] <= 0}}] Sp[t] (Ip[t] + Ipt[t] + Ic[t] + Ict[t]) +
            Piecewise[{{\[Iota]p, Mh[t] > 0}, {\[Iota]p2, 
               Mh[t] <= 0}}] ( h/(n h - h)) Sp[
             t] (Ip[t] + Ipt[t] + Ic[t] + Ict[t])) - \[Gamma]pt Ipt[t],
      Ict'[t] == 
       Piecewise[{{Tc, CT == 0 || Mh[t] > 0}, {0, 
            Mh[t] <= 
             0}}] (Piecewise[{{\[Beta]c, Mh[t] > 0}, {\[Beta]c2, 
               Mh[t] <= 0}}] Sc[
             t] (Ip[t] + Ipt[t] + Ic[t] + Ict[t]) + 
           Piecewise[{{\[Iota]c, Mh[t] > 0}, {\[Iota]c2, 
               Mh[t] <= 0}}] ( h/(n h - h)) Sc[
             t] (Ip[t] + Ipt[t] + Ic[t] + Ict[t])) - \[Gamma]ct Ict[t],
      
      IpC'[t] == 
       Piecewise[{{1 - Tp, CT == 0 || Mh[t] > 0}, {1, 
           Mh[t] <= 
            0}}] (Piecewise[{{\[Beta]p, Mh[t] > 0}, {\[Beta]p2, 
              Mh[t] <= 0}}] Sp[t] (Ip[t] + Ipt[t] + Ic[t] + Ict[t]) + 
          Piecewise[{{\[Iota]p, Mh[t] > 0}, {\[Iota]p2, 
              Mh[t] <= 0}}] ( h/(n h - h)) Sp[
            t] (Ip[t] + Ipt[t] + Ic[t] + Ict[t])),
      IcC'[t] == 
       Piecewise[{{1 - Tc, CT == 0 || Mh[t] > 0}, {1, 
           Mh[t] <= 
            0}}] (Piecewise[{{\[Beta]c, Mh[t] > 0}, {\[Beta]c2, 
              Mh[t] <= 0}}] Sc[t] (Ip[t] + Ipt[t] + Ic[t] + Ict[t]) + 
          Piecewise[{{\[Iota]c, Mh[t] > 0}, {\[Iota]c2, 
              Mh[t] <= 0}}] ( h/(n h - h)) Sc[
            t] (Ip[t] + Ipt[t] + Ic[t] + Ict[t])),
      IptC'[t] == 
       Piecewise[{{Tp, CT == 0 || Mh[t] > 0}, {0, 
           Mh[t] <= 
            0}}] (Piecewise[{{\[Beta]p, Mh[t] > 0}, {\[Beta]p2, 
              Mh[t] <= 0}}] Sp[t] (Ip[t] + Ipt[t] + Ic[t] + Ict[t]) + 
          Piecewise[{{\[Iota]p, Mh[t] > 0}, {\[Iota]p2, 
              Mh[t] <= 0}}] ( h/(n h - h)) Sp[
            t] (Ip[t] + Ipt[t] + Ic[t] + Ict[t])),
      IctC'[t] == 
       Piecewise[{{Tc, CT == 0 || Mh[t] > 0}, {0, 
           Mh[t] <= 
            0}}] (Piecewise[{{\[Beta]c, Mh[t] > 0}, {\[Beta]c2, 
              Mh[t] <= 0}}] Sc[t] (Ip[t] + Ipt[t] + Ic[t] + Ict[t]) + 
          Piecewise[{{\[Iota]c, Mh[t] > 0}, {\[Iota]c2, 
              Mh[t] <= 0}}] ( h/(n h - h)) Sc[
            t] (Ip[t] + Ipt[t] + Ic[t] + Ict[t])),
      
      Rp'[t] == \[Gamma]p Ip[t] + \[Gamma]pt Ipt[t],
      Rc'[t] == \[Gamma]c Ic[t] + \[Gamma]ct Ict[t],
      Mh'[
        t] == (PP (Sp[t] + Ipt[t] + Rp[t]) - 
          CC (Sc[t] + Ic[t] + Ip[t] + Ict[t] + Rc[t]) - 
          CT (Ipt[t] + Ict[t]))/n,
      Sp[0] == 0.4 n h - 1,
      Sc[0] == 0.6 n h,
      Ip[0] == 1,
      Ic[0] == 0,
      Ipt[0] == 0,
      Ict[0] == 0,
      IpC[0] == 0,
      IcC[0] == 0,
      IptC[0] == 0,
      IctC[0] == 0,
      Rp[0] == 0,
      Rc[0] == 0,
      Mh[0] == wealth
      },
     {Sp[t], Ip[t], Ipt[t], Rp[t], Sc[t], Ic[t], Ict[t], Rc[t], Mh[t],
       IpC[t], IptC[t], IcC[t], IctC[t]},
     {t, 0, tmax}
     ];
   Return[ODEDATA];
   ];