Use dimensionless variables $X=X_{\ast }x$, $Y=Y_{\ast }y$, $Z=Z_{\ast }z$, $t=T\tau$, withSubstitution gives the requested system, with dots now meaning differentiation with respect to $\tau$, andIn particular the same scale for $Y$ makes the first two coefficients of the cubic equal in magnitude.

Its derivative is $f^{\prime }(y)=1-2y+\gamma y^2$. For $0<\gamma <1$ its two simple zeros areThe cubic increases below $y_{-}$, decreases between the two fold points of a critical manifold, and increases above $y_{+}$. It has a local maximum at $y_{-}$ and a local minimum at $y_{+}$.

After a rapidly decaying transient response, $z\simeq f(y)$ when its relaxation is faster than the other variables. The reduced autonomous system isFor $ϵ\ll 1$, the fast $y$ dynamics at fixed $x$ attracts to the branches where $f^{\prime }(y)>0$ and repels from the middle branch. On the low branch $y<y_{-}<\alpha$, $x$ increases slowly. At its fold point of a critical manifold $x=f(y_{-})$, the attracting low branch disappears and $y$ jumps upward, with nearly constant $x$, to the high branch. On that branch $y>y_{+}>\alpha$, $x$ decreases slowly. At $x=f(y_{+})$ it jumps back to the low branch. These slow arcs and fast jumps form the singular relaxation cycle. For small positive $ϵ$, and parameters bounded away from the fold point of a critical manifold-endpoint cases, the nearby attracting periodic motion has the same pattern. The equilibrium $(f(\alpha ),\alpha )$ is on the repelling branch: its Jacobian matrix has determinant $1/ϵ>0$ and trace $-f^{\prime }(\alpha )/ϵ>0$.

The limiting landing values are $3/\gamma -2y_{-}$ after the upward jump and $3/\gamma -2y_{+}$ after the downward jump, obtained from the third root of $f(y)=f(y_{\pm })$. Thus the two vertical segments in the phase plane plot really connect roots of the same fast equation.

If $\gamma >1$, $f^{\prime }$ is strictly positive everywhere, so the cubic is monotone and has no fold points of a critical manifold. The reduced equilibrium is asymptotically stable. In factThe only invariant trajectory in the zero-derivative set is the equilibrium, so bounded reduced solutions converge to it; relaxation oscillations disappear.

The stronger condition $\boxed{\delta \ll ϵ}$ matters during a fast jump: $y$ then varies on the $ϵ$ time scale, not on the slow unit scale. Only if $z$ adjusts faster still may one keep $z\simeq f(y)$ through that jump. Merely imposing $\delta \ll 1$ would not justify this reduction on all pieces of the cycle. This is a three-time-scale cubic relaxation oscillator.

Write the equation as $ϵx^3-(x-1{)}^2=0$. For $ϵ\to 0^{+}$, two roots of a polynomial approach the double root one. Their splitting follows a fractional-power dominant balance, rather than an ordinary expansion in integer powers.

Put $\delta =\sqrt{ϵ}$ and $x=1+a\delta +b{\delta }^2+\cdots$. Substitution givesThus $a=\pm 1$ and $b=3/2$. The two bounded roots have the Puiseux series expansionsThe first two nonzero terms are $1\pm \sqrt{ϵ}$; the displayed correction also resolves their common order-$ϵ$ shift.

The third root escapes to infinity because the leading polynomial degree drops when $ϵ=0$. Balance $ϵx^3$ against $x^2$, giving $x=O({ϵ}^{-1})$. Set $X=ϵx$; the equation becomesFor $X=1+Aϵ+B{ϵ}^2+\cdots$, the first two coefficient equations are $A+2=0$ and $B+2A^2+2A-1=0$. HenceIts first two nonzero terms are ${ϵ}^{-1}-2$. These three distinct branches account for all roots of the cubic for sufficiently small positive $ϵ$. The Vieta formulas provide a consistency check: their sum is ${ϵ}^{-1}$, with the constant and order-$ϵ$ terms cancelling.

If $ϵ$ approaches zero through negative values, the bounded pair is complex: $x_{\pm }=1\pm i\sqrt{|ϵ|}+\frac{3}{2}ϵ+O(|ϵ{|}^{3/2})$. The large-root expansion is unchanged as a real asymptotic formula. The positive-parameter regime is used for the boundary-layer problem below.

Let $X=(x,y)$ be the integration variable and $P=(\xi ,\eta )$ the interior source point. Green's second identity, obtained from Green's theorem by applying the divergence theorem to $u\nabla G-G\nabla u$, isThe outward normal is used throughout. To obtain the prescribed representation with its positive volume term, chooseIndeed, substituting $\Delta u=f$ makes the left side $u(P)-{\int }_{D}Gf$. On the Dirichlet portion, $u=g$ and $G=0$ leave only $g{\partial }_{n}G$. On the Neumann portion, both normal derivatives vanish. Rearranging gives exactly the desired boundary-plus-volume representation. The distributional identity can equivalently be justified by excluding a small circle around $P$ and using the logarithmic singularity there.

For the quarter-plane, the fundamental solution of the Laplace equation with this sign convention is $\Phi (X,P)=(2\pi {)}^{-1}\log |X-P|$, since $\Delta \Phi ={\delta }_P$. The method of images uses an equal-sign reflection across the Neumann boundary $x=0$ and an opposite-sign reflection across the Dirichlet boundary $y=0$. HenceEquivalently,Only $(\xi ,\eta )$ lies inside the domain, so the interior Laplacian has exactly the required unit point source. At $y=0$ numerator and denominator agree, giving $G=0$. The formula is even in $x$, giving ${\partial }_{x}G=0$ at $x=0$, and therefore zero outward normal derivative there.

On the Dirichlet boundary $y=0$, the outward normal is $-{\partial }_y$, and differentiation gives the useful checkThis is the mixed Dirichlet-Neumann Green function in a quadrant. The image solution tends to zero at spatial infinity; on the unbounded domain the representation uses data and solutions for which the boundary term at infinity in Green's second identity vanishes.