Stability requirement for time step size for this discrete non-linear system

Question

The system of differential equations below evolve smoothly in time, however, when numerically evaluating the discretized system one must choose an appropriate time step. The required time step for stability depends on $\Delta$, the difference between the two boundary conditions: small values of $\Delta$ require small time steps. I would like to determine a relationship between the difference in the boundary conditions and the required time step.

Examining a linearized system indicates that the linearized system is stable for any size time step. This suggests that some non-linear analysis techniques are required. Any suggestions about where and how to analyze this system? Again the goal is to determine the maximum time step size $\Delta t$ which will allow the system to be integrated. (Note, this is not a homework problem.)

Below are the equations

$$ {dP_1\over dt} = k_1 \left(j_0 \, Sign\left(P_0 -P_{1(t)} \right)\sqrt{\left \vert P_0 -P_{1(t)}\right \vert} - j_1 \, Sign\left(P_{1 (t)} -P_{2(t)} \right)\sqrt{\left \vert P_{1(t)} -P_{2(t)} \right \vert }\right )$$

$$ {dP_2\over dt} = k_2\left(j_1 \, Sign\left(P_{1(t)} -P_{2(t)} \right)\ \sqrt{P_{1(t)} -P_{2(t)}} - j_2 \, Sign\left(P_{2(t)} -P_{3} \right)\sqrt{P_{2(t)} -P_{3}}\right )$$ with initial conditions $$P_0=P_3 + \Delta$$ $$P_{1(0)}=P_{2(0)}=P_3$$

The discretized system is below. $$P_1^{t+\Delta t} = P_1^{t} + \Delta t \, k_1\left(j_0 \, Sign\left(P_0 -P_1^t \right) \sqrt{ \left \vert P_0 -P_1^{t} \right \vert} - j_1 \, Sign\left(P_1^t -P_2^t \right) \sqrt{\left \vert P_1^{t} -P_2^{t} \right \vert }\right )$$

$$P_2^{t+\Delta t} = P_2^{t} + \Delta t \, k_2 \left(j_1 \, Sign\left(P_1^t -P_2^t \right) \sqrt{\left \vert P_1^t -P_2^{t} \right \vert} - j_2 \, Sign\left(P_2^t -P_3 \right)\sqrt{\left \vert P_2^{t} -P_3 \right \vert}\right )$$

---

Mathematica code

Below is some code to demonstrate the topic at hand. Choose a small enough time step and the results are OK. The question at hand is to estimate the required time step.

fStep[{ai_, aP1_, aP2_}] := Module[{bi, bP1, bP2, m0, m1, m2,
   \[CapitalDelta]t = 0.25 (1/2)^2, k1 = 2.5, k2 = 1, j0 = 23.3, 
   j1 = 23.3, j2 = 23.3, p0 = 152, p3 = 150 },
  (*note time step*)
  bi = ai + 1;
  (*compute flows given pressures*)
  m0 = Sign[p0 - aP1] j0 Sqrt[ Abs[p0 - aP1] ];
  m1 = Sign[aP1 - aP2] j1 Sqrt[ Abs[aP1 - aP2] ];
  m2 = Sign[aP2 - p3] j2 Sqrt[ Abs[aP2 - p3] ];
  (*compute pressures given flows*)
  bP1 = aP1 + \[CapitalDelta]t  k1 (m0 - m1);
  bP2 = aP2 + \[CapitalDelta]t k2 (m1 - m2);
  {bi, bP1, bP2}
  ]
r = NestList[fStep, {0, 150, 150}, 300];
Take[r, 10] // TableForm
Take[r, -10] // TableForm
vecP1 = Transpose[ {r[[All, 1]], r[[All, 2]]}];
vecP2 =  Transpose[ {r[[All, 1]], r[[All, 3]]}];
ListLinePlot[ {vecP1, vecP2}, PlotRange -> {{0, 40}, All}]