A characteristic curve in $(x,y,z)$ follows the vector field defining the first-order partial differential equation. Along such a curve the chain rule turns the equation intoIts characteristic projection is the curve $(x(s),y(s))$ in the independent-variable plane.

Writing $x(0)=\xi$, $y(0)=\eta$, $z(0)=\zeta$, separation of variables givesThese formulas hold on the maximal interval through $s=0$ where the denominators remain nonzero, using the real cube root. They include $\eta =0$ (a stationary characteristic projection) and $\zeta =0$ (the zero height). For $\zeta \ne 0$, the last formula follows from $d(z^{-3})/ds=1$. The alternative zero solution must be retained when deriving this by division by $z^4$.

At an equilibrium point, $x=-y^2$ and $x^2=a^{2}y$, so $y(y^3-a^2)=0$. With $d=|a{|}^{2/3}>0$, the two equilibrium points and their Jacobian matrix areAt $(-d^2,d)$, $\det J=-3a^2<0$. The trace-determinant stability criterion therefore gives an unstable saddle equilibrium. At the origin, the eigenvalues are $\lambda =\pm i|a|$: its linearization is a center equilibrium. The linearization stability theorem does not classify this nonhyperbolic equilibrium of the nonlinear system.

Nonlinear classification. In fact the answer depends on $|a|$:Here is an explicit derivation, rather than an inference from the imaginary eigenvalues. When $a^2=1$, set $L=1-x+y$. Near the origin $L>0$, and the first integralhasIts Hessian matrix at the origin is the identity. Small level curves are closed ovals on which the vector field is nonzero, proving the center equilibrium classification and Lyapunov stability.

For $a^2\ne 1$, a higher-order Lyapunov function calculation settles the issue. Put $k=|a|$, $U=x$, $V=ky$, $\tau =kt$, $p=k^{-1}$ and $q=k^{-2}$. ThenWrite $R^2=U^2+V^2$ and introduce the differential operatorsSuccessively cancel the terms of degrees three, four and five in the derivative of $R^2+H_3+H_4+H_5$, using these homogeneous polynomials:Direct differentiation verifies $\mathcal{L}{H}_3=-\mathcal{N}{R}^2$, $\mathcal{L}{H}_4=-\mathcal{N}{H}_3$ and $\mathcal{L}{H}_5=-\mathcal{N}{H}_4$. The remaining degree-six term has circular averageFor example, the only nonzero averages of degree-six monomials needed here are $⟨U^6⟩=⟨V^6⟩=5/16$ and $⟨U^4{V}^2⟩=⟨U^2{V}^4⟩=1/16$. Since $\mathcal{L}={\partial }_{\theta }$ in polar coordinates, the zero-average part can be cancelled by a degree-six homogeneous polynomial $H_6$ satisfyingIndeed its restriction to the unit circle is obtained by integrating a zero-average trigonometric polynomial; its Fourier modes have orders $2,4,6$, so multiplication by $R^6$ gives a polynomial. Consequently the positive Lyapunov function $W=R^2+H_3+H_4+H_5+H_6$ satisfiesFor $k>1$ this is strictly negative near the origin, giving asymptotic stability. For $k<1$ it is strictly positive, and $W_{\tau }\ge C{W}^3$ in a sufficiently small neighbourhood; every nonzero starting point must eventually leave that neighbourhood. Finally ${\theta }_{\tau }=1+O(R)>0$, so the nearby trajectories rotate: these are weak foci, with the stability stated above. Merely calling the nonlinear origin a center for every $a\ne 0$ would miss this distinction.

An equilibrium point $p=(x_0,y_0)$ of an autonomous planar system is a point where its vector field vanishes:For a locally well-posed ordinary differential equation, Lyapunov stability means that every prescribed neighbourhood can be protected by choosing a sufficiently small initial neighbourhood. More precisely, for every $\epsilon >0$ there exists $\delta >0$ such that every trajectory starting with $\Vert p(0)-p\Vert <\delta$ exists for all $t\ge 0$ and satisfiesThe number $\delta$ may depend on $\epsilon$ but not on the particular starting point. Asymptotic stability additionally requires convergence to the equilibrium point for all sufficiently close starting points.

Existence was proved by Picard iteration. If $y$ and $v$ are any two continuous solutions with the same $b\ne 0$, the previous part makes both bounded by $2B$. Subtracting their Volterra integral equations givesSince $\gamma <1$, this forces $y=v$.

For $b=0$, zero is a solution. On any finite interval $[0,T]$, another continuous solution has a finite bound $K_T={\max }_{[0,T]}|y|$. Since $K(t,s)\le 1$, its integral equation impliesThe Gronwall inequality gives $y=0$ on every such interval. Hence

Let $y$ be any continuous function solving the integral equation on $[0,\infty )$, and retain $B=|b|>0$. A first-exit argument avoids assuming in advance that its supremum norm is finite. If $|y|$ ever reached $2B$, continuity would give a first such time $T>0$. On $[0,T]$ we would have $|y|\le 2B$, so the kernel estimate would implya contradiction. Therefore $|y(t)|<2B$ for every finite $t$.

That pointwise strict inequality alone would not give a strict bound on the supremum norm. Apply the integral equation once more, now using $|y|\le 2B$ globally:

The difference of two successive Picard iterates factors asThe preceding bound and the kernel estimate therefore giveMoreover $D_1:=\Vert y_1-y_0{\Vert }_{\infty }\le B^2/C_1$. Summing the resulting geometric series, for $m>n$,Thus there is a uniform limit $y$ on $[0,\infty )$, continuous by the uniform limit theorem, with $\Vert y{\Vert }_{\infty }\le 2B$. The integral passes to this limit becauseTherefore $y$ satisfies the Volterra integral equation, and hence the initial value problem. This proves global existence by Global Picard iteration with a decaying kernel. The same estimate gives the useful error boundIf $b=0$, every iterate is zero and the convergence assertion is immediate.

Put $B=|b|>0$ and $\gamma =4B/C_1<1$. Initially $M_0=B<2B$. If $M_n<2B$, the preceding supremum norm estimate givesHence mathematical induction gives

For $s>t$, the coefficient bound gives $A(s)-A(t)\ge C_1(s-t)$. Thus for $0\le s\le t$ the positive kernel satisfiesThe starting Picard iterate is continuous and has supremum norm $M_0=|b|$. If the $n$th Picard iterate is bounded and continuous, its integral defines another continuous function, andTaking the supremum norm and applying mathematical induction proves boundedness of every iterate and

Write $A(t)={\int }_0^{t}a(s)ds$. Multiplying the ordinary differential equation by the integrating factor $e^{A(t)}$ and using the product rule givesIntegration from $0$ to $t$ therefore yields the Volterra integral equationConversely, a continuous function satisfying this integral equation becomes continuously differentiable after multiplication by $e^{A(t)}$, by the fundamental theorem of calculus, and satisfies the original initial value problem.

Fix any $t_0\ge 0$ and compare $y=R$ with $y=0$. The quotient required in a global Lipschitz condition in the dependent variable isConsequently,It is nevertheless locally Lipschitz in $y$: for $|y|,|v|\le R$,Thus the Picard-Lindelöf theorem gives local existence and uniqueness for the initial value problem; the remaining work establishes existence on the entire half-line for the specified small data.

One suitable autonomous planar system isIts equilibrium points are exactly the four combinations $(x,y)=(\pm 1,\pm 1)$. Its Jacobian matrix isThus $(-1,-1)$ and $(1,1)$ have eigenvalues $\pm 2$ and are saddle equilibria. At the other two points the eigenvalues are $\pm 2i$. To establish the nonlinear centre equilibrium classifications, use the first integralIts Hessian matrix is $\mathrm{diag}(2x,-2y)$, a negative-definite matrix at $(-1,1)$ and a positive-definite matrix at $(1,-1)$. Thus each is a strict local extremum, and nearby regular level sets are closed curves supporting periodic orbits, by the planar first integral criterion for a centre.

Immediately to the right of $(-1,1)$, the vector field points down because $x^2-1<0$, so that centre is clockwise. Immediately to the right of $(1,-1)$ it points up because $x^2-1>0$, so that centre is counterclockwise. This verifies both the nonlinear classifications and the prescribed orientations, with exactly four equilibrium points.

Applying time reversal of an autonomous differential equation to the system from part (i), and renaming its coordinates $(x,y)$, givesThe equilibrium points remain $(1,-1)$ and $(-1,-1)$. Time reversal of an autonomous differential equation preserves a saddle equilibrium, interchanging its stable and unstable directions; it also reverses the direction around a centre equilibrium. Therefore this system has precisely the required counterclockwise centre at $(1,-1)$ and saddle at $(-1,-1)$.

Using $x=-u$, $y=-v$ in the original autonomous planar system givesThe map $(x,y)\mapsto (-x,-y)$ is a rotation by $\pi$, preserving time and the clockwise or counterclockwise orientation of trajectories. Hence the given centre equilibrium at $(-1,1)$ becomes a clockwise centre at $(1,-1)$, and the given saddle equilibrium at $(1,1)$ becomes a saddle at $(-1,-1)$. These are the only equilibrium points, because rotation is invertible.

Set $z=y-\beta$. Then the autonomous planar system becomeswhich has no $\beta$ dependence. ThereforeMore precisely, replacing $\beta$ by $\beta +\delta$ sends every trajectory $(x(t),y(t))$ to $(x(t),y(t)+\delta )$, with the same time parametrization. The equilibrium points, nullclines, periodic orbits and separatrices all undergo this same translation. Their types and orientations are independent of $\beta$, so varying $\beta$ produces no bifurcation.

The nullclines areOn $y=\beta$, away from the equilibrium points, the vector field is vertical. Its crossing direction is downwards for $x<-1$, upwards for $-1<x<0$, downwards for $0<x<1$, and upwards for $x>1$, because $\dot{y}=x(x^2-1)$. On each of $x=-1,0,1$, the vector field is horizontal: it points right when $y>\beta$ and left when $y<\beta$.

Throughout each complementary region, $\dot{x}$ has the sign of $y-\beta$, and $\dot{y}$ has the four signs just listed. These signs orient the level sets from part (ii). In particular, the inner periodic orbits are clockwise. The saddle level is $h=1/4$, whose separatrices satisfyFor $-1<x<1$, the upper arc $y-\beta =(1-x^2)/\sqrt{2}$ goes from $(-1,\beta )$ to $(1,\beta )$; the lower arc goes back from $(1,\beta )$ to $(-1,\beta )$. These two heteroclinic orbits form a closed boundary around the family of periodic orbits. The remaining branches of the same saddle level extend outwards. The sketch shows the surrounding nonclosed trajectories as well as the nullclines and equilibrium points.

Where $y\ne \beta$, taking the quotient of the two ordinary differential equations givesChoose the first integralThe chain rule checks conservation even at points excluded in taking that quotient:Thus each trajectory is contained in a connected component of a level set of $h$. At $(0,\beta )$ the gradient vanishes and the Hessian matrix is the identity, so $h$ has a strict local minimum there. To make the closed-trajectory argument explicit, put $z=y-\beta$ and $0<E<1/4$. The bounded component of $h=E$ isThe two arcs join at $(\pm r_E,0)$ to form a simple closed curve that is a regular level set. The gradient cannot vanish on it: the only critical point of a dynamical system in $|x|<1$ is the centre, whose level is $0$.

The vector field is tangent and nowhere zero on this compact set, so its speed has a positive minimum. It therefore traverses the curve in finite time and returns to its initial point, giving a periodic orbit. These ovals shrink to the centre as $E\downarrow 0$. This is the planar first integral criterion for a centre and proves the required nonlinear assertion; it also gives Lyapunov stability without asymptotic stability.

The equilibrium points satisfy $y=\beta$ and $x(x^2-1)=0$, so they areThe Jacobian matrix and its eigenvalue equation areAt $x=\pm 1$ the eigenvalues are $\pm \sqrt{2}$, so both points are saddle equilibria. Their stable and unstable tangent directions have slopes $-\sqrt{2}$ and $\sqrt{2}$, respectively.

At $x=0$ the eigenvalues are $\pm i$, giving a centre for the linearized equation. Purely imaginary eigenvalues alone do not classify the nonlinear system. The first integral established in part (ii) proves that $(0,\beta )$ is in fact a centre equilibrium. It is clockwise: immediately to its right, with $0<x<1$ and $y=\beta$, the vector field points downwards. Thus the final nonlinear classification is

Continuous dependence of an ODE solution on parameters here means continuity of the map $\beta \mapsto y_{\beta }$ into $C([a-h,a+h])$ with its supremum norm: for every ${\beta }_0\in I$ and $\epsilon >0$, all sufficiently close $\beta \in I$ satisfySubtract the two integral equations and use the uniform Lipschitz constant $L$ on the rectangle. With $v(x)=|y_{\beta }(x)-y_{\gamma }(x)|$ we obtainPart (i), with $c=|\beta -\gamma |$ and $A=L$, therefore proves, for $L>0$,The case $\beta =\gamma$ follows from uniqueness or the same estimate with zero forcing. If $L=0$, the direct integral estimate gives $\Vert y_{\beta }-y_{\gamma }{\Vert }_{\infty }\le h|\beta -\gamma |$. In either case the map is even Lipschitz continuous, and hence continuous, throughout $I$ in the stated supremum norm.

Treat both sides of $a$ at once: for $\sigma \in \{-1,1\}$ set $w(t)=v(a+\sigma t)$, $0\le t\le h$. Nonnegativity givesDefine $W(t)=ct+A{\int }_0^{t}w(r)dr$. By the fundamental theorem of calculus, $W^{\prime }(t)=c+Aw(t)\le c+AW(t)$ and $W(0)=0$. The integrating factor $e^{-At}$ givesSince $w\le W$, this proves the affine Gronwall estimate

Near $(1,1)$ the solution of a differential equation has $y>0$, so the separable differential equation givesEquivalently its trajectory lies on the level setOn $y>0$, the function $G(y)=y-\log (1+y)$ is strictly increasing because $G^{\prime }(y)=y/(1+y)>0$, and its range is $(0,\infty )$. Thus the positive branch exists exactly whereThe connected component containing $x=1$ isThe implicit function theorem gives a smooth solution of a differential equation at every such $x$. It remains finite at every finite $x>x_{\ast }$, because $G$ increases without bound. As $x\downarrow x_{\ast }$, we have $y\downarrow 0$, and the original right side $x+x/y$ becomes singular. Hence this is the maximal interval $(x_{\ast },\infty )$ through the prescribed initial condition.

In particular, the solution is not defined for all $x\in (0,2)$. It cannot pass through $x=x_{\ast }$ as a classical solution, since its required value there is $y=0$.