(*Simulation Parameters and definitions*)
(*---------------------------------*)
Clear[f]
Clear[i, P, B]
f[P_, B_] := 1/2 P + 10 B/(1 + B);
tmax = 25;
eps = 0.001;
randNum = 100;
A = {{1/20, 1/4, 1/50}, {1/4, 1/26, 1/40}};
point = {16.666666666666735`, 0.`, 8.333333333333345`};

(*ODE System*)
ODEsys = {i'[t] == f[P[t], B[t]] - i[t],
   P'[t] ==
    P[t] (1 - A[[1, 1]] P[t] - A[[1, 2]] B[t] - A[[1, 3]] i[t]),
   B'[t] == B[t] (1 - A[[2, 2]] B[t] - A[[2, 3]] i[t])};

(*Generating parametric g curve.*)
(* Consider this section a  *)
(*--------------------------*)
set = List[];

While[Length[set] < randNum,
 holdSet =
  Join[set, Map[point + # &, RandomReal[{-eps, eps}, {randNum, 3}]]];
 set = Select[holdSet, #[[2]] >= 0 &];]

set = Drop[set, -(randNum - Length[set])];

(*Simulation*)
sol = ParametricNDSolveValue[{ODEsys, {P[0] == init1, B[0] == init2,
     i[0] == init0}}, {P, B, i}, {t, 0, tmax}, {init1, init2, init0}];

(*Averging solution over multiple inital conditions.*)
gCurve[t_] :=
  Evaluate[Mean[Through /@ (sol[#[[1]], #[[2]], #[[3]]][t] & /@ set)]];
gPlot = ParametricPlot3D[gCurve[t], {t, 0, tmax}, PlotRange -> All,
   ImageSize -> Large, PlotStyle -> Black];

(* Plot tubes of varying radiuses *)
(*----------------------------*)

(* Function which varies radius along the tube *)
radialfun[x_, radius_] := radius (Abs[x - 25]/50)^(1.20);

(* Controls the number of points plotted, 600+ gives a high quality
plot *)
numPoints = 50;

(* Controls opacity *)
a = 0.2;

{plot, vals} =
  Reap[ParametricPlot3D[gCurve[t], {t, 0, tmax},
    EvaluationMonitor :> Sow[t], MaxRecursion -> 0,
    Method -> {"TubePoints" -> numPoints}, PlotPoints -> numPoints,
    PlotStyle -> RGBColor[0., 0.5, 0.73, a] ]];
(* Plot different level sets of tubes*)
TubePlot1 =
  plot /. Line[pts_, rest___] :>
    Tube[pts, radialfun[vals[[1]], 1], rest];
TubePlot2 =
  plot /. Line[pts_, rest___] :>
    Tube[pts, radialfun[vals[[1]], 4], rest];
TubePlot3 =
  plot /. Line[pts_, rest___] :>
    Tube[pts, radialfun[vals[[1]], 10], rest];
TubePlot4 =
  plot /. Line[pts_, rest___] :>
    Tube[pts, radialfun[vals[[1]], 15], rest];

(*Show the plots of the different level tube sets*)
Show[gPlot, TubePlot1, TubePlot2, TubePlot3, TubePlot4,
 LabelStyle -> Automatic, BoxStyle -> Dashing[{0.02, 0.02}],
 Boxed -> None, PlotRangePadding -> 0, AxesOrigin -> Automatic,
 AxesEdge -> {{1, 1, 1}, None, Automatic}, Ticks -> Automatic,
 AxesStyle -> Thickness[0.005], ImageSize -> Large,
 ImageResolution -> Automatic, PlotRangeClipping -> True,
 PlotRange -> {{0, Automatic}, {0, Automatic}, {0, Automatic}}]