Replace backward Euler with BDF2

The time derivative is written literally as `(u - u_n) / dt` everywhere, so
FESTIM is first-order accurate in time. BDF2 is second order for the same cost
per step — one extra stored `Function`, same Jacobian sparsity, same Newton
solve — and much less diffusive on pulsed loading and stiff trapping
transients, where Euler damping currently forces small steps.

BDF2 is A- and L-stable like backward Euler, so this is a replacement, not an
opt-in scheme.

## Why it's small

Every term except the time derivative is already fully implicit at `t^{n+1}`
(reaction rates, `D(T)`, Sieverts/Henry BCs, surface reactions), so only the
derivative term changes. Same reason to prefer BDF2 over Crank-Nicolson, which
would need two-time-level assembly across `reaction.py`, the BCs and sources.

## What to change

**One shared helper** in `festim.helpers`, replacing five verbatim copies:

```python
form += time_derivative(u, u_n, u_nm1, dt, dt_prev) * v * dx(vol.id)
```

with, for `w = dt_n / dt_{n-1}`:

`du/dt = [ (1+2w)/(1+w) u^{n+1} - (1+w) u^n + w**2/(1+w) u^{n-1} ] / dt_n`

The variable-step form is mandatory — `dt` changes on adaptive stepping,
milestones and the final-time clip, so uniform coefficients would be wrong on
most real models.

**Call sites**: `HeatTransferProblem`
([:200](src/festim/heat_transfer_problem.py#L200)), `HydrogenTransportProblem`
([:935](src/festim/hydrogen_transport_problem.py#L935)), `...Discontinuous`
([:1625](src/festim/hydrogen_transport_problem.py#L1625)), `...ChangeVar`
([:2445](src/festim/hydrogen_transport_problem.py#L2445)), and the enclosure
pressure ODE ([:1927](src/festim/hydrogen_transport_problem.py#L1927)) — the
last must move with the rest or the coupled system drops back to first order.

**A `u_nm1` beside every `u_n`**: `HeatTransferProblem`
([:100](src/festim/heat_transfer_problem.py#L100)), `HydrogenTransportProblem`
([:732](src/festim/hydrogen_transport_problem.py#L732)) and
`spe.prev_solution`, `subdomain.u_n` +
`species.subdomain_to_prev_solution` ([:1538](src/festim/hydrogen_transport_problem.py#L1538)),
`GasSpecies.prev_solution`.

**A `dt_prev` constant** beside `self._dt`, and the shift in `iterate()`
([problem.py:248](src/festim/problem.py#L248),
[hydrogen_transport_problem.py:2339](src/festim/hydrogen_transport_problem.py#L2339))
done back-to-front. The `dt` history must shift *before* `modify_value`
overwrites `self.dt.value` — easy to get subtly wrong.

**First step** falls back to Euler (no `u^{n-1}` yet); one first-order step is
`O(dt^2)` globally, so second order is preserved. Internal detail, no user
switch, and initial conditions need no change. No restart needed after a
milestone or stepsize jump.

## Stepsize growth needs a study, not a textbook number

Variable-step BDF2 is zero-stable only for a *sustained* ratio
`dt_n/dt_{n-1} < 1+sqrt(2)`; isolated jumps are harmless.
`Stepsize.growth_factor` is unbounded ([stepsize.py:95](src/festim/stepsize.py#L95)),
so `growth_factor=3` would silently go unstable.

Before writing any guard: sweep `growth_factor` from 1.0 past 2.414 on a
closed-form case (exponential decay, or diffusion into a semi-infinite slab)
and plot error against realised step ratios. That gives what the bound can't —
where accuracy actually falls off versus where zero-stability is formally lost,
whether `max_stepsize`/milestone clipping already keeps realised ratios safe,
and whether failure is graceful or catastrophic (which decides warn vs clamp).
Note the configured `growth_factor` and the realised ratio are different checks
with different false-positive rates; the study should say which to guard.

## Risks of making it unconditional

Results change for every user. Not an API break, but needs a release note, and:

- **Positivity** — Euler is unconditionally monotone, BDF2 isn't, and can
undershoot negative near sharp fronts or BC step changes. Not cosmetic here:
negatives feed `sqrt` in the Sieverts law and can hard-fail Newton. Test a
step-change BC case before committing.
- **Tolerances** — system/MMS tests, `test/benchmark.py` and V-V-report cases
are all tuned against Euler. Expect a sweep of updates.

If undershoot is a real problem for trapping cases, fall back to a pluggable
scheme object on `Settings` (`F.BDF2()`, Euler default) — more machinery, wrong
default, but the helper above is where it would slot in.

## Testing

- **unit** — helper coefficients for uniform and non-uniform `dt` against
hand-computed values, Euler first step, the growth guard.
- **integration** — `u_nm1` created for every problem class; after two
`iterate()` calls the history and `dt` history hold the right values in order.
- **analytical** — the growth sweep above, kept as a test once the safe range
is known; also the step-change-BC positivity check.
- **system (MMS)** — halve `dt` on a fine mesh, assert observed order ~2, for
heat, hydrogen (continuous + discontinuous) and the enclosure ODE.

Run under `mpirun -n 2`; the shift is a local `x.array` copy, nothing special.

## Notes

- `CoupledHeatTransferHydrogenTransport` stays first order regardless — the
staggered heat->hydrogen ordering, not the scheme. Document, don't fix here.
- Follow-up issue: the BDF2/BDF1 difference is a nearly free local error
estimate, letting `Stepsize` adapt on accuracy instead of Newton iteration
count ([stepsize.py:144](src/festim/stepsize.py#L144)). Arguably the bigger win.

Start with the five derivative call sites in src/festim/heat_transfer_problem.py and src/festim/hydrogen_transport_problem.py, then inspect problem.py and stepsize.py for time and solution history updates. Review the unit, integration, analytical, and MMS tests described in the issue; done means variable-step BDF2 histories and first-step Euler behavior work across all listed systems, with the growth study, positivity check, and release-note impact addressed.

Written by the indexing model from the issue text.