Help debugging NDSolve with both a For loop and Table

Question

I am trying to run an ode system for a set number of time steps, take the size at the final time step for select equations and use these to calculate initial conditions for the next set of time steps. I am able to do this fine using a For loop, but now I would like to also vary one or more parameters. I have tried using both Table and Map to no avail. Any suggestions would be much appreciated. Here is a shortened version of my code:

lst = Table[X0 = 600;  
  Y0 = 300;  \[Mu]x = 0.006; \[Mu]y = 0.006; \[Mu]2 = 
   0.00001; \[Lambda] = 1; \[Sigma]x = 1;  \[Sigma]y = 1;  \[Gamma]x =
   0.1; \[Gamma]y = 0.2; \[Epsilon]y = 10; \[Phi] = 10; \[Alpha] = 
   0.008;
  gx[t_, \[Phi]_, \[Lambda]_, \[Epsilon]x_] := 
   If[t < \[Epsilon]x, 0, 
    PDF[GammaDistribution[\[Phi], 1/\[Lambda]], t - \[Epsilon]x]];
  gy[t_, \[Phi]_, \[Lambda]_, \[Epsilon]y_] := 
   If[t < \[Epsilon]y, 0, 
    PDF[GammaDistribution[\[Phi], 1/\[Lambda]], t - \[Epsilon]y]];
  soln[\[Epsilon]x_] =
   sol := NDSolve[{
      \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(x1[t]\)\) == 
       x3[t] gx[t, \[Phi], \[Lambda], \[Epsilon]x] - 
        x1[t] (\[Gamma]x + \[Mu]x),
      \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(x2[
         t]\)\) == \[Gamma]x x1[t] - \[Mu]2 x2[t],
      \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(y1[t]\)\) == 
       y3[t] gy[t, \[Phi], \[Lambda], \[Epsilon]y] - 
        y1[t] (\[Gamma]y + \[Mu]y),
      \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(y2[
         t]\)\) == \[Gamma]y y1[t] - \[Mu]2 y2[t],
      \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(x3[t]\)\) == 0,
      \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(y3[t]\)\) == 0,
      x1[0] == 0, x2[0] == 0, y1[0] == 0, y2[0] == 0, x3[0] == X0, 
      y3[0] == Y0}, 
     {x1, x2, x3, y1, y2, y3}, {t, 0, 120}];
  Xtotal = {X0}; Ytotal = {Y0};
  For[y = 0, y < 30, y++,
   loop = sol;
   X0 = \[Sigma]x y2[120] /. loop[[1, 4]];
   Y0 = \[Sigma]y x2[120]/(1 + \[Alpha] x2[120]) /. loop[[1, 3]];
   Xtotal = Append[Xtotal, X0];
   Ytotal = Append[Ytotal, Y0];
   ];, {\[Epsilon]x, 0, 10, 5}]
ListPlot[Xtotal]
ListPlot[Ytotal]

Answer 1

Taking into account our discussion, I propose the second variant depending on the parameters:

    tp = 120; \[Mu]x = 0.006; \[Mu]y = 0.006; \[Mu]2 = 
    0.00001; \[Lambda] = 1; \[Sigma]x = 1;  \[Sigma]y = 1;  \[Gamma]x =
    0.1; \[Gamma]y = 0.2; \[Epsilon]y = 10; \[Phi] = 10; \[Alpha] = 
    0.008;
  gx[t_, \[Phi]_, \[Lambda]_, \[Epsilon]x_] := 
     If[t < \[Epsilon]x, 0, 
       PDF[GammaDistribution[\[Phi], 1/\[Lambda]], t - \[Epsilon]x]];
  gy[t_, \[Phi]_, \[Lambda]_, \[Epsilon]y_] := 
     If[t < \[Epsilon]y, 0, 
       PDF[GammaDistribution[\[Phi], 1/\[Lambda]], t - \[Epsilon]y]];
     sol = ParametricNDSolve[{
          \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(x1[t]\)\) == 
            x3[t] gx[t, \[Phi], \[Lambda], \[Epsilon]x] - 
              x1[t] (\[Gamma]x + \[Mu]x),
          \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(x2[
       t]\)\) == \[Gamma]x x1[t] - \[Mu]2 x2[t],
          \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(y1[t]\)\) == 
            y3[t] gy[t, \[Phi], \[Lambda], \[Epsilon]y] - 
              y1[t] (\[Gamma]y + \[Mu]y),
          \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(y2[
       t]\)\) == \[Gamma]y y1[t] - \[Mu]2 y2[t],
          \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(x3[t]\)\) == 0,
          \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(y3[t]\)\) == 0,
          x1[0] == 0, x2[0] == 0, y1[0] == 0, y2[0] == 0, x3[0] == X0, 
          y3[0] == Y0}, 
        {x1, x2, x3, y1, y2, y3}, {t, 0, tp}, {\[Epsilon]x, X0, Y0}];
Table[{X0 = 600;  
     Y0 = 300;   Xtotal [\[Epsilon]x] = {X0}; 
   Ytotal [\[Epsilon]x] = {Y0};
     For[y = 0, y < 30, y++,
          {X0 = \[Sigma]x y2[\[Epsilon]x, X0, Y0][tp] /. sol;
        Y0 = \[Sigma]y x2[\[Epsilon]x, X0, Y0][
          tp]/(1 + \[Alpha] x2[\[Epsilon]x, X0, Y0][tp]) /. sol;
        Xtotal [\[Epsilon]x] = Append[Xtotal[\[Epsilon]x], X0];
        Ytotal [\[Epsilon]x] = Append[Ytotal[\[Epsilon]x], Y0];}
       ];}, {\[Epsilon]x, 5, 10, 5}];
Table[{ListPlot[Xtotal[\[Epsilon]x], PlotRange -> All, 
   PlotStyle -> Red, PlotLabel -> Row[{"\[Epsilon]x=", \[Epsilon]x}]],
   ListPlot[Ytotal[\[Epsilon]x], PlotRange -> All, 
   PlotStyle -> Orange]}, {\[Epsilon]x, 5, 10, 5}]