Implicit control via IFT/DAE — `Flow(ocp)` with `∂H̃/∂u = 0`

## Goal

The current `Flow(ocp, law)` constructor requires an explicit control feedback `law(x, p, v) -> u`. The goal is to support the case where the control is **not provided explicitly** but instead computed automatically from the pseudo-Hamiltonian stationarity condition `∂H̃/∂u = 0` — i.e., the smooth optimal control case derived from the PMP.

This is the missing branch of the with-control OCP flow story:

| Constructor | Control | Status |
|-------------|---------|--------|
| `Flow(ocp, law)` | explicit feedback | ✅ implemented |
| `Flow(ocp)` with implicit control | via `∂H̃/∂u = 0` | ❌ this issue |

Reference implementation in the old gallery: [gitlab.inria.fr/ct/gallery/examples/control-implicit](https://gitlab.inria.fr/ct/gallery/-/tree/master/examples/control-implicit). Roadmap: [v4 §1, Discussion #299](https://github.com/control-toolbox/CTFlows.jl/discussions/299).

---

## Two strategies

### Strategy A — IFT (Implicit Function Theorem)

Differentiate `∂H̃/∂u = 0` along the flow to get an explicit ODE for `u`:

```
u̇ = -(∂²H̃/∂u²)⁻¹ · (∂²H̃/∂z∂u)ᵀ · ẋ_H, z = (x, p), ẋ_H = (∂H̃/∂p, -∂H̃/∂x)
```

The augmented state `(x, p, u)` is integrated as an ODE. Requires an initial control guess `u(t0)`. The old-gallery RHS computes `∂H̃/∂z`, `∂H̃/∂u`, `∂²H̃/∂u²`, `∂²H̃/∂z∂u` by AD and sets:
- `dw[1:2n] = ẋ_H` (Hamiltonian vector field)
- `dw[2n+1:end] = -(∂²H̃/∂u²)⁻¹ · (∂²H̃/∂z∂u)ᵀ · ẋ_H`

**Pros:** pure ODE, no new heavy dependency. Prototype this first.

### Strategy B — DAE

Keep `∂H̃/∂u = 0` as an algebraic constraint and integrate `(x, p, u)` as a DAE (`DAEProblem`, IDA/Sundials or mass-matrix DFBDF).

**Pros:** no need for an initial `u(t0)` guess. **Cons:** new dependency (Sundials or similar), justifies its own weak-dep extension.

➡️ Prototype IFT first (pure ODE, no new dep), then decide whether DAE warrants the additional machinery.

---

## New AD primitives required in `CTBase.Differentiation`

Currently `CTBase.Differentiation` provides `pseudo_hamiltonian_gradient` / `pseudo_variable_gradient` (first-order). The IFT RHS needs:

- `∂²H̃/∂u²` — Hessian of `H̃` w.r.t. control `u`
- `∂²H̃/∂z∂u` — mixed Jacobian (derivative w.r.t. `z = (x,p)` of `∂H̃/∂u`)

These are new primitives to add alongside the existing ones. Default implementation: ForwardDiff on ForwardDiff (or ForwardDiff on the existing `PseudoHamiltonianSystem`).

---

## Call signature

For implicit-control flows, the initial control `u0` is a **keyword argument**, consistent with `variable` and `variable_costate`:

```julia
# point evaluation
xf, pf, uf = flow(t0, x0, p0, tf; control=u0, variable=v)

# trajectory
sol = flow((t0, tf), x0, p0; control=u0, variable=v)
```

`control` is **required** for implicit-control flows and **rejected** otherwise (trait-based validation, same pattern as `variable` for `NonFixed`/`Fixed`).

---

## Design decision: option vs. strategy

The implicit-control method can be surfaced as:
- A **keyword option** on `Flow(ocp)` / on the call (e.g. `implicit_control = :ift | :dae`), routed through the existing action-option machinery (like `hamiltonian_type`)
- A **`CTBase.Strategies` strategy** — first-class, dispatchable, pluggable (same family as `:di` backend and `:sciml` integrator strategies)

➡️ Decision pending. Lean toward a strategy if IFT/DAE needs its own options and extension boundary (DAE pulls in Sundials as a weak dep). **Decide together with #93** (flow derivative / IND) so both features use the same convention.

---

## Open questions

- IFT Hessian via AD (ForwardDiff) or user-supplied? Default: AD.
- Which DAE solvers to support, and behind which weak-dep extension?
- Interaction with `variable_costate` augmentation: does the augmented state become `(x, p, u, pv)` or `(x, p, pv, u)`? Order and coercion of the augmented block.

---

## Work items

- [ ] Decide option vs. strategy (together with #93)
- [ ] Prototype IFT RHS functor in `probe/` (augmented `(x, p, u)` ODE, validated against old gallery)
- [ ] Add `∂²H̃/∂u²` and `∂²H̃/∂z∂u` to `CTBase.Differentiation`
- [ ] Implement `Flow(ocp)` implicit-control path with `control=u0` keyword dispatch
- [ ] Trait-based validation: `control` required / rejected based on flow type
- [ ] Integration test: DI energy (`u* = p[2]` recovered implicitly from `∂H̃/∂u = 0`)
- [ ] (Optional) DAE path via `DAEProblem`, behind its own extension