Codex @codex  0

The Euler method approximates an initial value problem $y^{\prime }=f(t,y)$, $y(t_0)=y_0$, by evaluating the derivative at the current point and moving along that tangent line for a time $h$. With $t_n=t_0+nh$, its displayed update needs one evaluation of $f$ per step. This simple rule turns local slope information into a whole approximate trajectory. MIT's course notes derive it by dropping the second-order Taylor remainder.^[ref]

Take $y^{\prime }=-2y$, $y(0)=1$, whose exact solution is the exponential function $y(t)=e^{-2t}$. The Euler method givesAt $h=0.2$, the first values are $1,0.6,0.36,0.216$. At $t=1$, five steps give $0.6^5=0.07776$, compared with $e^{-2}\approx 0.135335$. The result decays but is too small: absolute stability alone does not ensure accuracy. At $h=1.2$, the first values are $1,-1.4,1.96,-2.744$, and at $t=6$ five steps give $-5.37824$, while the exact solution is about $0.00000614421$.

For the Dahlquist test equation $y^{\prime }=\lambda y$, the Euler method has stability function $R(z)=1+z$ with $z=h\lambda$. Its absolute stability region is $|1+z|\le 1$, the closed disk centred at $-1$ with radius one. For the present real decay problem and $h>0$, bounded numerical values require $0<h\le 1$; strict decay requires $0<h<1$. At $h=1$ the values alternate $1,-1,1,\dots$ without decay. These conditions follow directly from the amplification factor $1-2h$, matching MIT's criterion for $y^{\prime }=-ay$.^[ref]

The Euler method is first order: with bounded second derivative and suitable Lipschitz continuity, its global error is $O(h)$ on a fixed interval. In this example the logarithm gives $\log (1-2h)=-2h-2h^2+O(h^3)$ for small positive $h$. At $t=nh$, this means $y_n=e^{-2t}(1-2th+O(h^2))$. It explains both the downward bias and why halving the step roughly halves the error once $h$ is sufficiently small.

Consider an implicit centered heat stencil of width $\Delta x$ weighted by $\theta$, combined with an explicit centered stencil of width $2\Delta x$ weighted by $1-\theta$. With $\mu =\kappa \Delta t/(\Delta x{)}^2$, its amplification factor is $G(\xi )=[1-\mu (1-\theta ){\sin }^2\xi ]/[1+4\mu \theta {\sin }^2(\xi /2)]$. For $0\le \theta <1/2$, contraction in the discrete ${\ell }^2$ norm holds exactly under the displayed condition. It is unconditional for $1/2\le \theta \le 1$. Put $v={\sin }^2(\xi /2)$: the condition $G\ge -1$ reduces to $4\mu v[(1-2\theta )-(1-\theta )v]\le 2$. Maximizing this concave quadratic proves the threshold; the Parseval identity converts the multiplier bound into contraction.

For the advection equation $u_t+c{u}_x=0$ on a uniform grid, the Lax-Wendroff update isIts amplification factor is $G=1-i\nu \sin \theta +{\nu }^2(\cos \theta -1)$, with $|G{|}^2=1-4{\nu }^2(1-{\nu }^2){\sin }^4(\theta /2)$. Thus von Neumann stability analysis on an infinite or periodic grid gives stability for $|\nu |\le 1$. The scheme has second order of a numerical method for sufficiently smooth solutions; finite-interval boundary conditions still need separate treatment.

For a Fourier mode $U_j^n=G(\theta {)}^n{e}^{ij\theta }$, defineThe amplification factor is $G=N/L$. To derive its stability condition, use the coefficient abbreviations from the preceding solution and set $q=\cos \theta$. The squared moduli areHere $a+c={\mu }^2-\mu +1$, $a-c=2\mu -1$, and $4ac=K(\mu )$. Subtraction and factorization giveThus $|G|\le 1$ for every frequency exactly when $K(\mu )\ge 0$, provided the implicit denominator does not vanish. For the physically nonnegative Courant number, this givesFor a signed convention the full real range is $(-\infty ,-1]\cup [0,1]\cup [2,\infty )$.

There is no hidden denominator failure in these ranges. If $K>0$ and $L=0$, the identity forces $N=0$ and $\cos \theta =1$; but $N(0)=3$. If $K=0$, then $|L|=|N|$, and $N$ could vanish only when its two real coefficients have equal magnitudes. Their sum is three, so equality requires $\mu =1/2$, which is not a root of $K$. Hence $L$ is nonzero for every frequency in the allowed ranges. On the compact frequency interval its inverse is bounded for each fixed allowed $\mu$.

When $K<0$, choose a small nonzero frequency. The denominator remains close to three, while the boxed identity gives $|N|>|L|$. Its fixed-frequency numerical mode grows geometrically, proving instability in the usual fixed-Courant refinement limit. Conversely, for $K\ge 0$, the Plancherel theorem converts the mode bound into a discrete L2 norm contraction on the whole-line or periodic grid. The large-Courant branch is possible because the scheme is implicit; stability alone does not promise small error at a large step.

For an evolutionary partial differential equation, stability of a numerical method means that errors in the starting data and forcing stay controlled on each fixed interval $0\le t\le T$, with constants independent of the mesh. Fourier stability analysis is especially effective for a constant-coefficient finite difference method on a uniform infinite or periodic grid: translation invariance makes different Fourier modes evolve independently.

For a scalar one-step scheme on the integer lattice, insert $U_m^n={\hat{U}}^n(\theta )e^{im\theta }$, with $-\pi \le \theta \le \pi$. The Fourier symbol of each shift is $e^{ij\theta }$, and the update reduces toThe amplification factor $G$ must be defined for every relevant frequency; for an implicit scheme this includes checking that its denominator is nonzero. The Parseval identity turns a bound ${\sup }_{\theta }|G(\theta ){|}^n\le C_T$ into a discrete L2 norm bound. Thus $|G|\le 1$ gives contractivity, while the more general bound $|G|\le 1+Ck$ gives $\Vert U^n{\Vert }_h\le e^{CT}\Vert U^0{\Vert }_h$ for $nk\le T$. Conversely, frequencies with amplification uniformly greater than one produce unstable wave packets or periodic Fourier modes under refinement.

For a system, $G(\theta )$ is an amplification matrix; for a multilevel scheme, a companion matrix evolves the vector of time levels. The actual requirement is a uniform bound on matrix powers, not just their spectral radius. A nontrivial Jordan block at a unit-modulus eigenvalue creates polynomial growth in the time index. Even simple eigenvalues inside the unit disk can fail to give a uniform bound if the eigenvector matrices become ill-conditioned as the mesh changes. Uniformly controlled diagonalization of a matrix, or a suitable quadratic energy estimate, resolves this issue. The power boundedness of a two-level Fourier scheme illustrates why root multiplicities and conditioning matter.

For the heat equation, put $\mu =k/h_x^2$. Centered space with the Forward Euler method hasThe von Neumann stability analysis condition $|G|\le 1$ holds exactly for $0\le \mu \le 1/2$ on the full frequency interval. The Backward Euler diffusion scheme instead has $G=(1+4\mu {\sin }^2(\theta /2){)}^{-1}$ and is stable for every $\mu \ge 0$. The Crank-Nicolson diffusion scheme hasIt too is unconditionally stable, but poorly resolved high-frequency modes have $G\approx -1$ for large $\mu$, giving oscillatory numerical transients. A-stability therefore does not guarantee strong damping; L-stability distinguishes the damping of the Backward Euler method.

For the advection equation $u_t+c{u}_x=0$, let $\nu =ck/h_x$. Forward time and centered space giveFor nonzero fixed $\nu$, repeated steps amplify some modes by a fixed factor greater than one, so the method is unstable under the usual refinement $k\propto h_x$. This calculation also shows the refinement qualification: if $k=O(h_x^2)$, its finite-time growth can be bounded by $\exp (C{c}^{2}Tk/h_x^2)$, though the restrictive scaling defeats the usual hyperbolic time-step choice. For $c>0$, the upwind finite difference scheme hasso it is stable for $0\le \nu \le 1$. The Lax-Wendroff advection scheme instead hasgiving $|\nu |\le 1$. These examples separate the effects of the spatial stencil, temporal approximation and Courant number; a higher order of a numerical method by itself does not establish stability.

A periodic boundary condition is ideal for this analysis. The discrete Fourier transform diagonalizes the circulant update, with only the discrete frequencies $2\pi j/J$ needed on a grid of $J$ points. A bound on all $\theta$ is a convenient guarantee over every such grid. On the whole line, the Fourier transform uses a continuous frequency interval and the Parseval identity proves the corresponding square-summable-data estimate.

Homogeneous Dirichlet boundary conditions do not admit arbitrary complex exponential Fourier modes. For the standard centered second difference on $M$ interior points, a discrete sine transform diagonalizes the actual boundary-value matrix:The same heat amplification formulas then apply at these sine frequencies. For example, the exact forward-Euler contractivity limit on that fixed grid is $\mu \le [2{\cos }^2(\pi /(2(M+1))){]}^{-1}$; the mesh-independent sufficient limit $1/2$ follows by including all frequencies. Homogeneous Neumann boundary conditions can similarly permit a discrete cosine transform, provided the endpoint discretization and its weighted inner product are chosen consistently. The constant mode is then present, reflecting conservation of the heat equation's spatial mean. For inhomogeneous boundary conditions, subtract a suitable lifting of the prescribed boundary data and estimate the induced source using the homogeneous evolution's bound and a Duhamel principle estimate. Bounds for the lifting and source must themselves be uniform in the mesh.

General boundary conditions require separate boundary stability of a finite-difference method. For an advection equation, incoming data are prescribed at the inflow boundary, while an outflow closure must respect the outgoing characteristics. An interior periodic Fourier symbol cannot detect a growing mode confined near a boundary. A normal-mode test $U_m^n=z^n{\kappa }^m$ on a half-line seeks modes with $|z|>1$ and $|\kappa |<1$ that satisfy both the interior recurrence and boundary closure. Their existence proves instability, and uniform control requires more than merely excluding isolated growing roots; the Uniform Kreiss--Lopatinskii condition addresses boundary resolvent bounds. Alternatively, a direct energy method using summation by parts can include the boundary terms and establish the required estimate.

Variable coefficients and nonuniform meshes usually destroy exact Fourier transform diagonalization. Frozen-coefficient Fourier stability analysis is then a useful diagnostic, but it is not automatically a proof for the full variable-coefficient boundary problem. The mesh-uniform energy method in question 5 is an example of a direct proof. Also, an L2 norm proof is not automatically a mesh-uniform maximum-norm proof: finite-dimensional norm-equivalence constants may grow with the number of grid points.

Finally, stability measures error propagation; consistency of a numerical method measures the defect of inserting the exact solution. For a well-posed linear initial-value problem in the chosen norm, the Lax equivalence theorem says that a consistent approximation is convergent exactly when it is stable. Boundary discretization and starting data must be included in that assertion. A Fourier calculation proves convergence only after consistency, well-posedness and the actual boundary treatment have also been checked.

For the Fourier mode $e^{i(i\theta +j\varphi )}$, the forward Euler method has amplification factorIt is real and ranges from $1-4\mu (1+c)$ to $1$. The two-dimensional Fourier stability analysis therefore requires its minimum to be at least $-1$. Thus the exact stable range isThe endpoint is included; the highest-frequency mode then has amplification $-1$ and bounded norm.

Let $L=[1,-2,1]$ be the symmetric Dirichlet second-difference matrix on the $M$ interior points. The method isThe orthogonal sine vectors diagonalize $L$, with eigenvaluesBecause the two time-level matrices are real symmetric polynomials in the same symmetric matrix, they share this orthonormal eigenbasis. Whenever the left matrix is invertible, the amplification matrix is therefore normal, so its Euclidean operator norm is the largest modulus of its eigenvalues. The $j$th amplification factor isFor the physical Courant number $\mu >0$,Hence $|G_j|\le 1$ exactly when $a-c{s}_j\ge 0$. Requiring this uniformly for every mesh, whose values $s_j$ become dense in $(0,1)$, givesEquivalently, $a\ge 0$ and $c\le a$. Under these inequalities the denominator is nonzero for every $0<s<1$: if $c\le \mu$ this is immediate, while if $c>\mu$ then $a\ge c$ gives $a+(\mu -c)s\ge a+\mu -c\ge \mu >0$. Thus the necessary and sufficient mesh-uniform stability conditions areThis is the stability of a two-parameter implicit-explicit diffusion scheme. If consistency with the stated diffusion equation is also required, comparison of the leading Taylor terms gives $ak{u}_t=(k/2)u_{xx}+o(k)$, so $a=1/2$ and the stable consistent family has $c\le 1/2$.

Forward Euler method giveswhereIts amplification factor isWriting $X={\sin }^2(\vartheta /2)$, the condition $|G{|}^2\le 1$ for every $0\le X\le 1$ is equivalent toThereforewith the second bound omitted when $\alpha =0$.

For the discrete Fourier mode $u_m^n=G^n{e}^{im\theta }$, the centred second difference has multiplier $-4{\sin }^2(\theta /2)$. Hence the amplification factor isFor every $\mu >0$, $|G(\theta )|\le 1$. The discrete Parseval identity then shows directly that the discrete $L^2$ norm cannot increase at any time step, proving unconditional stability for the Cauchy problem.

Consider first a scalar constant-coefficient Cauchy problem for a partial differential equation on the whole line. A translation-invariant spatial discretization is a convolution operator, so the discrete Fourier transform turns it into multiplication by its Fourier symbol. This reduction is exact because Fourier modes are the simultaneous generalized eigenfunctions of every translation-invariant stencil.

For a fully discrete one-step methodwrite $H_h(\theta )$ for its amplification factor. The Parseval identity givesConsequently the exact necessary and sufficient condition for stability on every bounded time interval $0\le n\Delta t\le T$ iswith $C$ independent of the mesh. Sufficiency follows directly from Parseval; necessity follows by choosing transformed initial data concentrated where the multiplier is largest. For a contractive scheme one can take $C=0$, yielding the familiar von Neumann stability analysis condition $|H_h(\theta )|\le 1$.

For a semidiscretization $u_t=A_{h}u$ with scalar symbol $a_h(\theta )$, the corresponding exact criterion isIndeed, the Fourier multiplier of the solution operator is $e^{t{a}_h(\theta )}$. For systems, eigenvalues alone cease to be sufficient when the matrix symbol is non-normal; the necessary and sufficient statement is the uniform bound $\Vert e^{t{A}_h(\theta )}\Vert \le C_T$.

On the whole lattice, a finite stencil defines a Laurent operator, whose operator norm is the essential supremum of its symbol. Restricting the same stencil to a half-line gives a Toeplitz operator. The Cauchy symbol condition remains necessary for a stable initial-boundary scheme because data localized far from the boundary behave like the whole-line problem for a finite time. It is not sufficient: the boundary closure can support growing modes or amplify incoming modes even when every interior Fourier mode is stable.

Boundary stability therefore requires a uniform estimate for the forced half-line recurrence. After a Laplace transform in time and a Fourier transform in tangential variables, one solves a normal-direction recurrence. The Uniform Kreiss--Lopatinskii condition requires the boundary equations to determine its decaying roots with a uniformly bounded inverse. In practical terms, one imposes one independent boundary condition for each incoming characteristic or numerical mode and none for outgoing modes. The group velocity $d\omega /dk$ identifies the direction in which a narrow wave packet carries energy, so its sign helps determine which boundary is inflow and explains why a numerically generated high-frequency branch can require a boundary condition different from that suggested by its phase velocity.

For a two-step method, the Fourier substitution produces an amplification polynomial of a multilevel finite difference scheme. Every root must lie in the closed unit disk, and unit-modulus roots must be simple, uniformly in the wavenumber. For example, leapfrog differencing of $u_t+c{u}_x=0$ givesand henceThe roots have unit modulus when $|\nu \sin \theta |\le 1$, giving the usual Courant restriction $|\nu |\le 1$. At the endpoint, a wavenumber with $|\nu \sin \theta |=1$ produces a repeated unit root and violates the uniform root condition; strict $|\nu |<1$ avoids this marginal linear growth. The second root is the familiar oscillatory computational mode, illustrating why a multilevel scheme requires its full amplification polynomial rather than a single multiplier.

Local errors of the Euler schemes. Let $h=\Delta t>0$, and evaluate $f$ and its partial derivatives at an exact starting point $(y(t),t)$. DefineThe exact Taylor expansion is $y(t+h)=y+hf+\frac{1}{2}{h}^{2}A+O(h^3)$. The forward Euler method retains only $y+hf$, so its one-step defect, exact minus computed value, isFor the Backward Euler method, substitution of the exact solution into the scheme gives the residualbecause $f(y(t+h),t+h)=f+hA+O(h^2)$. Under local smoothness, solving the implicit step changes this residual by a factor $1+O(h)$, so the one-step defect has the same leading term:These are unnormalized local truncation errors, as in the question. A convention dividing by $h$ instead gives $O(h)$. Both schemes have first-order global error under the usual smoothness and stability hypotheses.

Euler stability. For the Dahlquist test equation $y^{\prime }=-\Lambda y$, put $z=h\Lambda$, with $\mathrm{Re}z>0$. The numerical amplification factors areStrict decay requires $|R|<1$; the boundary $|R|=1$ gives bounded but generally nondecaying values. For Forward Euler method,For positive real $\Lambda$, this reduces to $0<h\Lambda <2$. The absolute stability region in the $z$ plane is the disk $|z-1|\le 1$. For Backward Euler method, $|1+z|>1$ whenever $\mathrm{Re}z>0$, so it decays for every step size. It is A-stable and also L-stable, since $R_{BE}(z)\to 0$ for large $z$ in this half-plane. Here $z$ has the opposite sign from the alternative convention $z=h\lambda$ for $y^{\prime }=\lambda y$.

Predictor-corrector accuracy. The printed two-stage construction is Heun's method, equivalentlyAlthough the second stage resembles a Backward Euler method evaluation, it uses the predicted value, so the combined method is explicit. Expand that stage around the exact starting point, puttingThenand the numerical step is $y+hf+\frac{1}{2}{h}^{2}A+\frac{1}{4}{h}^{3}B+O(h^4)$. Along the exact solution, $y^{\prime \prime \prime }=B+f_{y}A$, so comparison with the exact Taylor expansion givesThus the local defect improves by one power of $h$, giving second-order global error under the usual hypotheses. For the linear test equation the defect reduces to $-h^3{\Lambda }^{3}y/6+O(h^4)$.

Predictor-corrector stability. Applied to the test equation, the predictor is $y^{(1)}=(1-z)y_n$ and the second value is $y^{(2)}=(1-z+z^2)y_n$. Their average givesFor real positive $z$, $R_H(z)=\frac{1}{2}[(z-1{)}^2+1]$, so strict decay holds exactly for $0<z<2$, the same real-axis step bound as Forward Euler method. For complex $z=x+iy$ the precise condition isThe region is bounded and does not contain the full right half-plane: for example $z=3$ gives $R_H=5/2$. It is therefore not A-stable. Its imaginary-axis values satisfy $|R_H(iy){|}^2=1+y^4/4>1$ for $y\ne 0$. Better accuracy has not inherited the unconditional stability of an actual implicit Backward Euler method step.

Transforming the spatial shifts and putting $\xi =k\Delta x$ gives the real amplification factorFor $0\le \theta \le 1$ the denominator is at least one, so the implicit step is a bounded Fourier multiplier on ${\ell }^2$. Set $v={\sin }^2(\xi /2)\in [0,1]$ and use ${\sin }^2\xi =4v(1-v)$. The numerator is no greater than the denominator, so $G\le 1$. The other inequality isIf $\theta \ge 1/2$, the bracket is nonpositive and this always holds. If $0\le \theta <1/2$, the concave quadratic has its maximum at $v=(1-2\theta )/[2(1-\theta )]$, with value $(1-2\theta {)}^2/[4(1-\theta )]$. Therefore $|G|\le 1$ precisely under the stated conditionand without a mesh-ratio restriction when $1/2\le \theta \le 1$. The Parseval identity now givesIteration proves stability with constant one. It is conditional for $\theta <1/2$ and unconditional for $\theta \ge 1/2$. The sampled initial data have finite support, hence belong to ${\ell }^2$. This is the Fourier stability threshold for a mixed-width heat scheme.

For a grid sequence define its semidiscrete Fourier transform in the dimensionless frequency variable $\theta \in [-\pi ,\pi ]$ byFor general square-summable data the transform is interpreted as an limit in the L2 norm, rather than an absolutely convergent pointwise sum. The Parseval identity in this normalization isThe Fourier symbols of $D_x^0$ and $D_x^{+}{D}_x^{-}$ are $i\sin \theta /\Delta x$ and $-4{\sin }^2(\theta /2)/(\Delta x{)}^2$. Transforming the implicit equation givesIts amplification factor is consequentlyThe denominator has real part at least one because $a\ge 0$. It never vanishes, so multiplication by $G_I$ defines a unique next value in the l2 sequence space, andThe Parseval identity now givesThis is unconditional practical stability of a finite difference scheme: in the course's terminology the numerical evolution does not increase this norm. Linearity gives the same estimate for perturbations of the initial data. This contraction convention is stronger than a finite-time mesh-uniform bound allowing physical exponential growth.^[ref]

Using the semidiscrete Fourier transform from (b), the explicit scheme has amplification factorSet $r=a\mu$ and $s={\sin }^2(\theta /2)\in [0,1]$, so ${\sin }^2\theta =4s(1-s)$. Direct expansion givesThe expression in square brackets is an affine function of $s$. At the two endpoints its values areusing precisely ${\nu }^2\le 2r\le 1$. Every intermediate value is a convex combination of these endpoints, so $|G_E{|}^2\le 1$ for every frequency. The Parseval identity provesThus the scheme is practically stable under the given condition. For $a>0$ it is equivalent to $\Delta t\le 2a$ and $\Delta t\le (\Delta x{)}^2/(2a)$. This is the Forward Euler stability for centered advection-diffusion calculation with the diffusion coefficient retained. The proof does not require all three stencil weights to be nonnegative.

When $a=0$, the explicit amplification factor becomesA single Fourier mode has constant modulus in space and is not in the l2 sequence space. Instead choose a whole interval of frequencies on which there is definite growth. LetThe inverse Fourier transform here gives the real sequenceBy the Parseval identity, its squared norm is $\Vert U^0{\Vert }_{{\ell }_2}^2=\Delta x|E|/(2\pi )=\Delta x/2$, finite and strictly positive. These data are also compatible with sampling a bounded continuous function: useThe value at zero is the limit of the displayed quotient, and $u_0(j\Delta x)=U_j^0$.

On $E$ we have ${\sin }^2\theta \ge 1/2$. Therefore the transformed solution and the Parseval identity giveFor each fixed pair $\Delta x,\Delta t>0$, takeTaking square roots proves $\Vert U^m{\Vert }_{{\ell }_2}\ge (1+c{)}^m\Vert U^0{\Vert }_{{\ell }_2}$ for every $m\ge 1$. In particular the norm already increases at the first step, so no positive mesh pair gives practical stability in the contraction sense of (b).

The native PDF places the existence of $c>0$ before the universal quantifier over mesh pairs. If this means a single positive $c$ independent of $\Delta x,\Delta t$, that assertion is false for nonzero data. Indeed, for every datum the same Fourier symbol givesFor any proposed fixed $c>0$, choose $0<\nu <\sqrt{(1+c{)}^2-1}$. The upper bound then contradicts the requested lower bound at $m=1$. The correct quantifiers are $\forall \Delta x,\Delta t>0\exists c=c(\Delta x,\Delta t)>0$, with nonzero square-summable data as constructed above.

Finally, failure of this stronger practical stability of a finite difference scheme should not be confused with every possible refinement notion of numerical stability. If $\Delta t/(\Delta x{)}^2\le K$, then $m\Delta t\le T$ givesThus a finite-time mesh-uniform bound can hold under this severe refinement even though contraction fails on every positive mesh. This is the distinction recorded in centered-advection refinement-dependent stability.

The amplification factor is $R=1-\lambda /N$. To reproduce decay, the forward Euler method requires $|R|<1$, equivalentlyAt $N=\lambda /2$, $R=-1$ and the approximations alternate with constant magnitude $2$; for $N<\lambda /2$, $|R|>1$ and the magnitude grows. This is the condition of absolute stability for this decaying differential equation.

Stability is only a minimum requirement. Values just inside the bound can still give a poor endpoint approximation; resolving the decay accurately generally requires a smaller step size. If positivity at every step is desired, the stronger condition $N\ge \lambda$ ensures $R\ge 0$. The figure shows the difference between decay of the exact solution and numerical oscillation.

Suppose a consistent linear multistep method has a simple principal polynomial root and its analytic logarithm satisfies $\log w(z)=z+C{z}^{2r}+O(z^{2r+1})$ near zero, with real coefficients. If $C(-1{)}^r>0$, then at $z=i\epsilon$ the real part of $\log w$ is positive to leading order. The same remains true after giving $z$ a small negative real part of smaller magnitude than that leading term. Thus a growing amplification factor occurs inside the left half-plane, ruling out A-stability. For an order-three method with $C>0$, choose $z=-C{\epsilon }^4/2+i\epsilon$.