Boundary value problems — shooting method, finite difference: Worked examples

Before we start, a self-contained toolkit — every symbol used below, defined right here.

Definition The linear-BVP finite-difference template (defined here, used everywhere below)

A general linear second-order BVP is written $y^{''} = p (x) y^{'} + q (x) y + r (x), y (a) = α, y (b) = β .$

== $p (x)$ == multiplies the first derivative $y^{'}$ — the convection term (a "drift" that pushes one direction).
== $q (x)$ == multiplies $y$ itself — the reaction/stiffness term.
== $r (x)$ == is the source — everything with no $y$ in it (forcing, heating).

Chop $[a, b]$ into $N$ equal steps of width $h = \frac{b - a}{N}$ ; interior nodes are $x_{1}, \dots, x_{N - 1}$ . Replacing $y^{'}$ and $y^{''}$ by central differences and multiplying by $h^{2}$ gives, at each interior node $i$ , the three-coefficient stencil: $lower (1 + \frac{h}{2} p_{i}) y_{i - 1} diag - (2 + h^{2} q_{i}) y_{i} + upper (1 - \frac{h}{2} p_{i}) y_{i + 1} = RHS h^{2} r_{i} .$ where $p_{i} = p (x_{i})$ , etc. This is the only FDM formula we use; each example just plugs numbers into it.

Two more words we lean on constantly:

Interior node: a grid point strictly between the two ends — the ones we must solve for. The ends are already given.
Residual: how far our current answer is from satisfying the rule. For shooting it is $ϕ (s) = y (b; s) - β$ : the gap between where our shot landed and the bucket $β$ .

The scenario matrix

Cell	What makes it different	Example
A. Linear + shooting	$ϕ (s)$ is a straight line → exact in 2 shots	Ex 1
B. Nonlinear + shooting	$ϕ (s)$ is curved → must iterate secant	Ex 2
C. Linear FDM, $p = 0$	symmetric tridiagonal, no convection	Ex 3
D. Linear FDM, $p \neq = 0$	asymmetric off-diagonals (convection)	Ex 4
E. Degenerate grid ( $N = 2$ )	only ONE interior node → single equation	inside Ex 3
F. Limiting case: refine $h$	check error falls $\propto h^{2}$ (halve $h$ → $\frac{1}{4}$ error)	Ex 5
G. Sign check: which side did we land?	overshoot vs undershoot, negative residual	Ex 6
H. Real-world word problem	heated rod, units, physical read-off	Ex 7
I. Exam twist: derivative (Neumann) boundary	$y^{'} (b)$ given, not $y (b)$ → ghost node	Ex 8

Every method-cell (A–I) is hit at least once below.

Ex 1 — Cell A: Linear shooting, exact in two shots

Worked example Solve

$y^{''} = 4 y, y (0) = 0, y (1) = sinh 2$ by shooting. True solution is $y = sinh (2 x)$ , so $y^{'} (0) = 2$ . Pretend we don't know that.

Forecast: the ODE is linear, so how many shots do we need before the answer is exact — 2, 5, or "infinitely many"?

Shot 1: guess $s_{0} = 0$ . Solving $y^{''} = 4 y, y (0) = 0, y^{'} (0) = 0$ gives $y \equiv 0$ , so $y (1; 0) = 0$ . Why this step? We need one data point $(s, landing)$ . Zero slope is the cheapest probe — it tells us the baseline landing.
Shot 2: guess $s_{1} = 1$ . By linearity $y = \frac{1}{2} sinh (2 x)$ (slope-1 solution), so $y (1; 1) = \frac{1}{2} sinh 2 \approx 1.8134$ . Why this step? A line needs two points. Now we know $ϕ$ at two slopes.
Interpolate to the root. Target $β = sinh 2 \approx 3.6269$ : $s^{⋆} = s_{0} + (s_{1} - s_{0}) \frac{β - y ( 1 ; s _{0} )}{y ( 1 ; s _{1} ) - y ( 1 ; s _{0} )} = 0 + 1 \cdot \frac{3.6269 - 0}{1.8134 - 0} = 2.0.$ Why this step? For linear ODEs $y (b; s)$ is a straight line in $s$ , so the exact root is one interpolation away — no iteration.

Verify: $s^{⋆} = 2$ matches the true $y^{'} (0) = 2$ exactly. Forecast answer: 2 shots. ✓

Ex 2 — Cell B: Nonlinear shooting, secant must iterate

Worked example Solve

$y^{''} = y^{2}, y (0) = 0, y (1) = 1$ by shooting. Because of the $y^{2}$ the ODE is nonlinear, so $ϕ (s)$ is a curve, not a line.

Forecast: will two shots nail it, or will the residual still be nonzero after interpolating twice?

For a clean, checkable arithmetic we use the secant formula on made-up-but-consistent landing data. Suppose our IVP solver returns: $y (1; s_{0} = 0.5) = 0.60, y (1; s_{1} = 1.0) = 1.30.$ Target $β = 1$ , so residuals $ϕ_{0} = 0.60 - 1 = - 0.40$ , $ϕ_{1} = 1.30 - 1 = + 0.30$ .

Secant step. $s_{2} = s_{1} - ϕ_{1} \frac{s _{1} - s _{0}}{ϕ _{1} - ϕ _{0}} = 1.0 - 0.30 \cdot \frac{1.0 - 0.5}{0.30 - ( - 0.40 )} .$ Why this step? We approximate the curve $ϕ (s)$ by the straight line through the last two points and hit its root — the secant idea, no derivative of $ϕ$ needed.
Arithmetic. $\frac{0.5}{0.70} = 0.7143$ ; $0.30 \times 0.7143 = 0.2143$ ; so $s_{2} = 1.0 - 0.2143 = 0.7857$ . Why this step? This is the next, better slope guess.
Not done! We must fire a new IVP at $s_{2}$ , get a new residual, and repeat — because the true $ϕ$ curves away from our line. Why this step? On a curved $ϕ$ , one secant line only approximates the root; iteration drives $∣ ϕ ∣ \to$ tol.

Verify: the residual straddled zero ( $ϕ_{0} < 0 < ϕ_{1}$ ), so the root is bracketed between $0.5$ and $1.0$ ; our $s_{2} = 0.7857$ lands inside — a sane secant step. Forecast answer: iterate, not two shots. ✓

Ex 3 — Cells C & E: Linear FDM, $p = 0$ , degenerate single node

Worked example

$y^{''} = y, y (0) = 0, y (1) = sinh 1$ , take $N = 2$ so $h = 0.5$ . Match to the template: $p = 0, q = 1, r = 0$ . With $N = 2$ there is exactly one interior node $x_{1} = 0.5$ — the degenerate case.

Forecast: with only one unknown, how many equations do we solve — a system, or a single line?

Write the grouped stencil using the template with $p = 0, q = 1$ : $(1) y_{0} - (2 + h^{2}) y_{1} + (1) y_{2} = 0.$ Why this step? Substituting $p = 0$ kills the off-diagonal asymmetry (lower $=$ upper $= 1$ ); $q = 1, h^{2} = 0.25$ fixes the middle coefficient.
Insert boundary values $y_{0} = 0, y_{2} = sinh 1 = 1.1752$ : $0 - 2.25 y_{1} + 1.1752 = 0.$ Why this step? $y_{0}, y_{N}$ are known numbers — move them to the RHS (Trap "RHS").
Solve the one equation. $y_{1} = \frac{1.1752}{2.25} = 0.5223.$ Why this step? One interior node ⇒ the tridiagonal system collapses to a single scalar equation (Cell E).

Verify: true $sinh (0.5) = 0.5211$ ; error $\approx 0.0012$ . Small even with 2 intervals, and it will shrink like $h^{2}$ (proved next). ✓

Ex 4 — Cell D: Linear FDM with convection $p \neq = 0$

Worked example

$y^{''} = 2 y^{'} - y + e^{x}, y (0) = 1, y (1) = e,$ take $N = 4$ so $h = 0.25$ . Write the equation at node $i = 1$ ( $x_{1} = 0.25$ ). Match to the template: $p = 2, q = - 1, r = e^{x}$ . Now the off-diagonals differ — that asymmetry is the convection.

Forecast: will the coefficient multiplying $y_{i - 1}$ (lower) be bigger or smaller than the one on $y_{i + 1}$ (upper)?

Compute the three coefficients with $h = 0.25$ :
- lower $= 1 + \frac{h}{2} p = 1 + 0.25 = 1.25$ ,
- diag $= - (2 + h^{2} q) = - (2 + 0.0625 \cdot (- 1)) = - (2 - 0.0625) = - 1.9375$ ,
- upper $= 1 - \frac{h}{2} p = 1 - 0.25 = 0.75$ . Why this step? Direct substitution into the template's lower/diag/upper slots.
RHS $= h^{2} e^{x_{1}} = 0.0625 e^{0.25} = 0.0625 \cdot 1.2840 = 0.08025$ . Why this step? $r (x) = e^{x}$ evaluated at the node, scaled by $h^{2}$ .
Equation at $i = 1$ (with $y_{0} = 1$ moved right): $- 1.9375 y_{1} + 0.75 y_{2} = 0.08025 - 1.25 (1) = - 1.16975.$ Why this step? $y_{0}$ is known; subtract $1.25 y_{0}$ to the RHS (Trap "RHS").

Verify: with $p = 2 > 0$ , the term $\frac{h}{2} p$ adds to the lower ( $1.25$ ) but subtracts from the upper ( $0.75$ ), so lower $>$ upper — the upwind side is heavier, exactly what positive convection produces. Forecast answer: lower is bigger. ✓

Worked example Take the FDM of Ex 3 (

$y^{''} = y$ , exact $sinh$ ) at $x = 0.5$ for $N = 2$ and $N = 4$ . Confirm the error quarters.

Forecast: halving $h$ (from $0.5$ to $0.25$ ) should shrink the error by a factor of about... 2, 4, or 8?

$N = 2$ error (Ex 3): $∣0.5223 - 0.5211∣ = 0.00122$ . Why this step? This is our coarse-grid benchmark.
Set up the $N = 4$ system explicitly. Now $h = 0.25, h^{2} = 0.0625$ ; interior nodes $x_{1} = 0.25, x_{2} = 0.5, x_{3} = 0.75$ . With $p = 0, q = 1$ every diagonal is $- (2 + h^{2}) = - 2.0625$ and every off-diagonal is $1$ . Boundaries $y_{0} = 0, y_{4} = sinh 1 = 1.1752$ go to the RHS (Trap "RHS"), so $⎩ ⎨ ⎧ - 2.0625 y_{1} + y_{2} = - y_{0} = 0, y_{1} - 2.0625 y_{2} + y_{3} = 0, y_{2} - 2.0625 y_{3} = - y_{4} = - 1.1752.$ Why this step? A finer grid means three coupled interior equations — a genuine tridiagonal system.
Solve by the Thomas method (forward-eliminate, back-substitute).
- Row 1: $y_{1} = y_{2} /2.0625$ .
- Substitute into row 2: $(y_{2} /2.0625) - 2.0625 y_{2} + y_{3} = 0 \Rightarrow - 1.5776 y_{2} + y_{3} = 0 \Rightarrow y_{3} = 1.5776 y_{2}$ .
- Substitute into row 3: $y_{2} - 2.0625 (1.5776 y_{2}) = - 1.1752 \Rightarrow y_{2} (1 - 3.2538) = - 1.1752 \Rightarrow y_{2} = \frac{- 1.1752}{- 2.2538} = 0.52140.$ Why this step? Thomas is just Gaussian elimination that exploits the tridiagonal shape — eliminate below the diagonal, then back-substitute for the value we want, $y_{2} = y (0.5)$ .
$N = 4$ error: $∣0.52140 - 0.52110∣ = 0.00030$ . Ratio $\frac{0.00122}{0.00030} \approx 4.1$ . Why this step? This is the observed error-shrink factor.

Verify: the central stencil is $O (h^{2})$ , so halving $h$ predicts $\times \frac{1}{4}$ error. We saw $\approx 4.1 \times$ smaller — matches (Trap "refine": always refine to trust). Forecast answer: factor 4. ✓

Ex 6 — Cell G: Sign of the residual (over/undershoot)

Figure — Boundary value problems — shooting method, finite difference

Figure: three shooting attempts for $y^{''} = y, y (0) = 0$ , all launched from the left end. The magenta curve (slope $s_{0} = 0$ ) stays flat and ends below the navy dashed "bucket" line $β = sinh 1$ — an undershoot. The violet curve (slope $s_{1} = 2$ ) rises steeply and ends above the bucket — an overshoot. The orange curve (slope $s^{⋆} = 1$ ) lands exactly on the bucket. Read the sign of the residual straight off the picture: below the dashed line ⇒ $ϕ < 0$ , above ⇒ $ϕ > 0$ , on it ⇒ $ϕ = 0$ .

Worked example For

$y^{''} = y, y (0) = 0, y (1) = sinh 1 \approx 1.1752$ , interpret the two shooting shots by sign.

Forecast: if we shoot with slope $s_{0} = 0$ , do we land short of or beyond the bucket $β$ ?

Shot $s_{0} = 0$ : $y (1; 0) = 0$ . Residual $ϕ_{0} = 0 - 1.1752 = - 1.1752 < 0$ . Why this step? A negative residual means we landed below/short of the target — undershoot. In the figure this is the magenta path ending below the navy dashed bucket line.
Shot $s_{1} = 2$ : $y (1; 2) = 2 sinh 1 \approx 2.3504$ . Residual $ϕ_{1} = 2.3504 - 1.1752 = + 1.1752 > 0$ . Why this step? A positive residual means we shot over — the violet path ends above the bucket.
Root is bracketed. Since $ϕ_{0} < 0 < ϕ_{1}$ , the true slope lies between $0$ and $2$ ; interpolation lands at $s^{⋆} = 1$ (the orange path, which touches the bucket). Why this step? Opposite signs guarantee a crossing — the geometric meaning of "aim between a short shot and a long shot."

Verify: $s^{⋆} = 0 + 2 \cdot \frac{1.1752 - 0}{2.3504 - 0} = 1.0$ , the true slope. Sign logic and root agree. ✓

Ex 7 — Cell H: Real-world word problem (heated rod)

Worked example A rod of length

$L = 1 m$ conducts heat with an internal source. Steady temperature $T (x)$ obeys $T^{''} = - 100 (^{\circ} C / m^{2}), T (0) = 2 0^{\circ} C, T (1) = 2 0^{\circ} C .$ Estimate the mid-rod temperature $T (0.5)$ by FDM with $N = 2$ ( $h = 0.5$ ).

Forecast: with both ends at $2 0^{\circ}$ and heating inside, is the middle hotter or cooler than $2 0^{\circ}$ ?

Match to the template. Rewrite $T^{''} = 0 \cdot T^{'} + 0 \cdot T + (- 100)$ , so $p = 0, q = 0, r = - 100$ . Why this step? Reading off $p, q, r$ lets us reuse the one stencil formula.
Stencil at $x_{1} = 0.5$ , $h^{2} = 0.25$ : $y_{0} - 2 y_{1} + y_{2} = h^{2} r = 0.25 (- 100) = - 25.$ Why this step? $p = q = 0$ gives the pure "1, −2, 1" stencil; the RHS carries the source $r$ .
Insert $T_{0} = T_{2} = 20$ : $20 - 2 T_{1} + 20 = - 25 \Rightarrow - 2 T_{1} = - 65 \Rightarrow T_{1} = 32. 5^{\circ} C .$ Why this step? Both ends known; single equation gives the midpoint (Trap "RHS" applied).

Verify: exact solution $T = 20 + 50 x (1 - x)$ gives $T (0.5) = 20 + 50 (0.25) = 32. 5^{\circ} C$ — matches exactly (for constant $r$ the central stencil is exact on this parabola). Units are $^{\circ}$ C throughout. Forecast answer: hotter, $32. 5^{\circ}$ . ✓

Ex 8 — Cell I: Exam twist, a derivative (Neumann) boundary

Worked example

$y^{''} = y, y (0) = 0, y^{'} (1) = cosh 1 \approx 1.5431$ (a slope given at the right end). Take $N = 2$ , $h = 0.5$ ; set up the equations. The right boundary is now a derivative, not a value — the classic exam curveball.

Forecast: at $x = 1$ we have no $y_{2}$ value handed to us. Where does the missing equation come from?

Introduce a ghost node $y_{3}$ at $x_{3} = 1.5$ (just outside the rod) and write the central derivative at the boundary node $x_{2} = 1$ : $y^{'} (1) \approx \frac{y _{3} - y _{1}}{2 h} = cosh 1 \Rightarrow y_{3} = y_{1} + 2 h cosh 1 = y_{1} + 1.5431.$ Why this step? Using a central (not one-sided) difference at the boundary keeps $O (h^{2})$ accuracy (Trap "order"). The ghost node $y_{3}$ is the fictitious value the central formula needs.
Write the $y^{''} = y$ stencil at BOTH unknown nodes $i = 1, 2$ (here $p = 0, q = 1, h^{2} = 0.25$ , diag $= - 2.25$ ): $i = 1 : y_{0} - 2.25 y_{1} + y_{2} = 0, i = 2 : y_{1} - 2.25 y_{2} + y_{3} = 0.$ Why this step? Node 2 is now an unknown (only its slope was given, not its value), so it needs its own stencil equation — and that equation contains the ghost $y_{3}$ .
Eliminate the ghost. Substitute $y_{3} = y_{1} + 1.5431$ (from step 1) into the node-2 equation: $y_{1} - 2.25 y_{2} + (y_{1} + 1.5431) = 0 \Rightarrow 2 y_{1} - 2.25 y_{2} = - 1.5431.$ Why this step? The ghost value is now expressed in terms of a real unknown plus known data, so it disappears and we are left with a solvable $2 \times 2$ system.
Solve the $2 \times 2$ with $y_{0} = 0$ : from $i = 1$ : $- 2.25 y_{1} + y_{2} = 0 \Rightarrow y_{2} = 2.25 y_{1}$ . Sub into step-3 equation: $2 y_{1} - 2.25 (2.25 y_{1}) = - 1.5431 \Rightarrow (2 - 5.0625) y_{1} = - 1.5431 \Rightarrow y_{1} = \frac{- 1.5431}{- 3.0625} = 0.5039.$ Then $y_{2} = 1.1338.$ Why this step? Back-substitution — the Thomas idea by hand for a tiny system.

Verify: exact $y = sinh x$ gives $sinh (0.5) = 0.5211$ and $sinh (1) = 1.1752$ . Our coarse $N = 2$ estimates $0.5039, 1.1338$ are in the right ballpark and would converge as $h \to 0$ . The ghost-node trick correctly supplied the missing equation. ✓

Recall Rapid self-test

Which matrix cell needs iteration even though shooting is used? ::: Cell B — nonlinear ODE makes $ϕ (s)$ curved. A negative residual $ϕ (s_{0}) < 0$ means we shot… ::: short (undershot the target $β$ ). Halving $h$ in a central-difference FDM changes the error by roughly… ::: a factor of $\frac{1}{4}$ (it is $O (h^{2})$ ). With a Neumann boundary $y^{'} (b)$ given, what do we add? ::: a ghost node plus a central-difference equation for the slope. Why is the lower off-diagonal bigger than the upper when $p > 0$ ? ::: the convection term $\frac{h}{2} p$ adds to the lower, subtracts from the upper.

Worked examples — Boundary value problems — shooting method, finite difference

The scenario matrix

Ex 1 — Cell A: Linear shooting, exact in two shots

Ex 2 — Cell B: Nonlinear shooting, secant must iterate

Ex 3 — Cells C & E: Linear FDM, $p = 0$ , degenerate single node

Ex 4 — Cell D: Linear FDM with convection $p \neq = 0$

Ex 5 — Cell F: Limiting refinement, error $\propto h^{2}$

Ex 6 — Cell G: Sign of the residual (over/undershoot)

Ex 7 — Cell H: Real-world word problem (heated rod)

Ex 8 — Cell I: Exam twist, a derivative (Neumann) boundary