Shooting Method and bvp4c

1. The problem type

A typical reduced system: f‴ + ff″ − f′² − Mf′ = 0 with f(0) = 0, f′(0) = 1 and f′(∞) = 0, possibly coupled to θ″ + Pr fθ′ = 0 with θ(0) = 1, θ(∞) = 0. Conditions at two ends make it a boundary-value problem (BVP); the infinite domain is truncated to [0, η∞] with η∞ chosen large enough that the answers stop changing.

2. Shooting: turn the BVP into an IVP

Write the ODE as a first-order system y₁ = f, y₂ = f′, y₃ = f″. The initial condition y₃(0) = f″(0) = s is unknown — so guess it, integrate to η∞ with Runge–Kutta, and measure the miss:

Now the whole problem is root-finding on F(s): update s by the secant method or Newton's method and repeat until |F| is below tolerance. For the MHD problem the converged s must equal −√(1+M) — the exact solution grades the code.

3. Collocation: bvp4c

MATLAB's bvp4c instead represents the solution by piecewise cubics on a mesh and enforces the ODE at collocation points, solving the resulting nonlinear algebraic system globally with adaptive mesh refinement. It is more robust than shooting when the system is stiff, when parameters push solutions toward blow-up, or when dual solutions exist — the second branch is captured by supplying a different initial guess function, which is how stability-analysis papers locate both branches.

4. Validation discipline

Three checks appear in every sound paper: (i) reproduce a known exact value — e.g. f″(0) = −1 at M = 0, or Blasius's 0.33206; (ii) show insensitivity to η∞ and mesh refinement; (iii) tabulate against previously published values in limiting cases. A result that skips these is a plot, not a computation. The same standard drives this site's rule that every calculator engine is tested against known-correct values.

5. A genuine numerical experiment

Setup. Solve f‴ + ff″ − f′² − Mf′ = 0, f(0) = 0, f′(0) = 1, f′(η∞) = 0 by RK4 (step 0.001) plus root-finding on s = f″(0). Exact answers for grading: s = −1 (M = 0), s = −2 (M = 3). What follows are actual computed numbers, not illustrations.

Run 1 — it works. M = 3, η∞ = 10, secant iteration from s = −1.5, −2.5: converged to f″(0) = −2.000000, matching the exact −√(1+M) to six decimals. When shooting behaves, it is fast and clean.

Run 2 — the trap. M = 0 with a short domain η∞ = 6: bisection settles on s = −1.128 — 13% wrong, with a tiny residual proudly displayed. The culprit: demanding f′ = 0 exactly at finite η∞, when the true solution still carries f′(6) = e⁻⁶ ≈ 0.0025, forces the solver to a spurious root where f′ crosses zero on its way to divergence. At η∞ = 4 the same procedure gives −1.004 (small domain bias only); pushed to η∞ = 10 naively, the bracket collapses entirely.

Run 3 — why large η∞ breaks single shooting. Linearise about the solution: perturbations grow like e^{+βη}, so an error ε in s reaches the far boundary amplified by e^{βη∞} — at β = 1, η∞ = 10 that is e¹⁰ ≈ 22 000. The residual function F(s) becomes a cliff: astronomically steep near the root, spurious zeros elsewhere. This is textbook ill-conditioning of single shooting on problems with growing modes, demonstrated live.

The morals (each is standard professional practice): (i) target |f′(η∞)| < tol with tol ~ e^{−βη∞}, not exact zero; (ii) repeat at several η∞ and report insensitivity; (iii) when growth defeats you, switch to multiple shooting (subdivide the domain, match at interfaces — growth amplification per segment shrinks to e^{βΔη}) or to collocation. This experiment is the reason bvp4c dominates this literature.

6. bvp4c in practice

A complete, honest MATLAB skeleton for the M-problem:

with boundary residuals [y₁(0), y₂(0)−1, y₂(η∞)] and an initial guess like y = [η(1−e^{−η}) approximants]. Calls: solinit = bvpinit(linspace(0,10,50), @guess); sol = bvp4c(@ode, @bc, solinit); then read f″(0) = sol.y(3,1). Collocation solves globally — no growing-mode amplification — and the residual control adapts the mesh where the solution bends. Dual solutions (shrinking sheets, opposing buoyancy) are found by feeding a second, qualitatively different guess: the solver converges to whichever branch's basin the guess sits in, which is exactly how stability-analysis papers harvest both.

7. Common misconceptions

"Small residual = correct answer." Run 2's wrong answer had a beautiful residual. Residuals certify that an equation was solved on that domain — validation against exact limits and η∞-sweeps certify it was the right problem.

"η∞ should be as large as possible." For single shooting, larger η∞ can be strictly worse (Run 3). For collocation it is merely wasteful. The correct value is "large enough that answers stop moving," demonstrated, not asserted.

"bvp4c is a black box that can't be wrong." It faithfully solves what you pose: a wrong guess finds the other branch; too-loose tolerances smear f″(0) in the third decimal; a careless mesh near steep Hartmann-type layers misses them. The validation trio of §4 applies to collocation exactly as to shooting.

"Shooting is obsolete." For well-conditioned problems it is the fastest method there is, trivially parallel over parameters, and self-graded when exact limits exist — this site's own numerics used it. Know its failure mode (you have now seen it) and it remains a fine tool.

8. Practice problems

P1. Recast θ″ + Pr fθ′ = 0 as a first-order system for solution alongside the flow. Solution: y₄ = θ, y₅ = θ′; y₄′ = y₅, y₅′ = −Pr y₁y₅; BCs y₄(0) = 1, y₄(η∞) = 0 — a fourth and fifth equation appended to the momentum trio.

P2. With β = 2 (M = 3), how small is the neglected tail f′(η∞) at η∞ = 6? At 10? Solution: e^{−12} ≈ 6×10⁻⁶; e^{−20} ≈ 2×10⁻⁹ — set shooting tolerances accordingly, not tighter.

P3. Estimate the error-amplification factor of single shooting for M = 8 at η∞ = 10, and the per-segment factor if the domain is split into 5 for multiple shooting. Solution: β = 3: e³⁰ ≈ 10¹³ (hopeless); per segment e⁶ ≈ 403 (manageable) — multiple shooting's raison d'être in one line.

P4. Design the validation table for a new Casson–MHD paper. Solution: rows: (β_c → ∞, M = 0) vs Crane's −1; (β_c → ∞, M = 3) vs −2; grid halving; η∞ ∈ {8, 10, 12}; and comparison against one published table in the overlap case — the referee's checklist, pre-empted.

9. Going deeper

Modern practice adds continuation: solve at an easy parameter value, then march the parameter in small steps, feeding each solution as the next guess — the standard route to fold points where dual solutions merge. At the folds themselves, pseudo-arclength continuation parametrises along the solution curve and walks around the turning point that defeats naive marching; eigenvalue analysis of the linearised operator then labels the branches stable/unstable, completing the physical story the twin solutions began. These techniques — implemented in AUTO, MATCONT, or a page of your own MATLAB — are the professional finish on the numerics this chapter opened.

10. Historical context

"Shooting" is honest artillery language: guess the elevation, fire, observe the miss, correct — ranging a gun and solving a BVP are the same feedback loop, and the numerical method inherited the name early in the twentieth century. Its systematic theory (including multiple shooting, developed for aerospace trajectory work) was consolidated in Keller's classic monograph (1968) and in Stoer & Bulirsch. Collocation grew alongside; the specific tool that conquered this literature, MATLAB's bvp4c, was written by Kierzenka and Shampine (2001) — a three-stage Lobatto IIIa collocation scheme with residual-based mesh adaptation, chosen deliberately to be forgiving of the rough initial guesses users actually supply. Its siblings (bvp5c; SciPy's solve_bvp, a direct descendant) run the same design. The method's cultural dominance in stretching-sheet research is thus partly technical merit and partly the standard-tool effect: once benchmark tables were bvp4c tables, every successor validated against them with the same instrument.

11. Another way to see it: Newton shooting with sensitivities

Secant iteration treats the residual F(s) = f′(η∞; s) as a black box. Newton's method opens it: the derivative dF/ds can be computed exactly by integrating the variational (sensitivity) system alongside the flow. Differentiate the ODE with respect to the unknown initial slope s and write ξ(η) = ∂f/∂s: it obeys the linearised equation

Integrate the six-equation system (three for f, three for ξ) once per iteration; then F′(s) = ξ′(η∞) and Newton updates s ← s − F(s)/F′(s) with quadratic convergence — typically three or four integrations to machine precision from a sensible guess. The same construction generalises to systems (a Jacobian of sensitivities) and is exactly the machinery inside multiple-shooting codes. Conceptually it reframes shooting as root-finding with an exact derivative oracle: the ODE itself tells you how the target moves when the trigger moves — the artillery metaphor upgraded with a ballistics computer.

12. Frequently asked questions

How do I pick the initial guess for bvp4c? Feed it the physics you already know: the M = 0 exact solution (1 − e^{−η}, e^{−η}, −e^{−η}) is a superb guess for moderate M; for extended models, continuation from a solved neighbour beats cleverness every time.

Stiff versus ill-conditioned — the difference? Stiffness is an integrator's problem: widely separated decay rates force tiny steps on explicit methods (cure: implicit solvers). Shooting's disease is different — a growing mode amplifying initial-condition error (cure: multiple shooting or collocation). This chapter's Run 3 is the second disease; boundary-layer ODEs are rarely stiff in the first sense.

RK4 fixed-step or adaptive RK45? Adaptive for production (error control per step); fixed-step for the transparent experiments of §5, where you want the method's behaviour undisguised. Either way, halve the tolerance/step once and demand the answer stand still — the grid-independence ritual.

My two bvp4c runs give two different answers. Bug? Possibly a blessing: check whether both satisfy the equations (residuals) — if so you have found dual solutions, and the parameter regime (shrinking? opposing buoyancy?) will usually confirm you are near a fold. Report both; analyse stability; that is a paper section, not an error message.

13. Further practice

P5. Verify the sensitivity equation of §11 by differentiating f‴ + ff″ − f′² − Mf′ = 0 term by term with respect to s. Solution: ∂/∂s of each term with ξ = ∂f/∂s gives ξ‴ + ξf″ + fξ″ − 2f′ξ′ − Mξ′ = 0 ✓, and the initial data follow from differentiating f(0) = 0, f′(0) = 1, f″(0) = s.

P6. For the exact M = 0 solution, solve the sensitivity system's large-η behaviour and show F′(s)|_{s=−1} ~ growth like e^{+η∞}. Solution: the linearised operator about f = 1 − e^{−η} has far-field modes e^{−η}, e^{−η}·poly, and e^{+…η} from the f-coefficient asymptotics; the growing mode dominates ξ′(η∞) — Newton's derivative is huge, steps are tiny, and the cliff-shaped residual of Run 3 is re-derived analytically.

P7. Design a continuation path to reach a hybrid-nanofluid, M = 5, suction S = 1 solution from scratch. Solution: start (M, S, φ) = (0,0,0) with Crane's exact profile → march M: 0→5 in steps of 1 → march S: 0→1 in 0.25s → switch on φ to targets. Each stage seeds the next; any step that fails is halved. This recipe, verbatim, is how the tables in your thesis-adjacent papers are actually produced.

14. Worked exam problem

Problem. Set up — completely, ready to code — the coupled flow–heat problem for M = 2, Pr = 7: equations, first-order system, boundary conditions, initial guess, and the validation row you will publish.

Solution. Equations: f‴ + ff″ − f′² − Mf′ = 0; θ″ + Pr fθ′ = 0. System: y₁ = f, y₂ = f′, y₃ = f″, y₄ = θ, y₅ = θ′ with y′ = (y₂, y₃, y₁y₃−... sign care: y₃′ = −y₁y₃ + y₂² + My₂; y₄′ = y₅; y₅′ = −Pr y₁y₅). Boundary conditions: y₁(0) = 0, y₂(0) = 1, y₄(0) = 1; y₂(η∞) = 0, y₄(η∞) = 0, η∞ = 10. Guess: the exact flow at this M — β = √3, f = (1−e^{−βη})/β — plus θ = e^{−2η}. Validation row: the code must return f″(0) = −√3 = −1.7321 to four decimals (exact), and −θ′(0) should land near the √Pr estimate 0.798×√7 ≈ 2.11 (same order, somewhat below it at moderate Pr). A solver that passes both is entitled to its new-physics columns; one that cannot is not — this half-page is the methods section of a publishable paper.

15. Key takeaways

Shooting converts a BVP to iterated IVPs and is excellent when well-conditioned — but growing modes amplify guess errors by e^{βη∞}, spawning spurious roots and razor-thin brackets: this chapter demonstrated the failure with real numbers. Collocation (bvp4c) solves globally and dominates the literature for cause. Either way, validation is non-negotiable: exact limits, grid refinement, η∞ insensitivity — residuals certify solving, not solving the right problem.

16. Where to go next

The exact solutions used as graders come from the similarity chapter; the thermal equation being appended is derived in heat transfer; the coefficients for nanofluid runs come from nanofluid models. The Reynolds number calculator supplies instant exact values for any M — a pocket validation oracle.