Randomized stable SDIRK sampler

Integrate this at the trajectory-step operation, not inside the score or vector-field network. For a state z_n at time t_n, draw x once per step from the coefficient distribution used by the paper; for an MVP, use x uniformly on [0.25,0.5] and sweep this range to verify algebraic stability. Compute k1 by solving k1=f_theta(t_n+xh,z_n+hx k1). Compute k2 by solving k2=f_theta(t_n+(1-x)h,z_n+h[(1-2x)k1+xk2]), then return z_{n+1}=z_n+h(k1+k2)/2. Use two or three Newton iterations with Jacobian-vector products, or Anderson-accelerated fixed-point iteration; do not form a dense Jacobian. A Newton linear solve should use conjugate gradients, GMRES, or a low-rank preconditioner for the operator I-hx J_theta, matching the paper's resolvent structure. For diffusion models, apply the method to the reverse-time probability-flow ODE; for neural ODE training, use it in the forward solver with checkpointing or differentiation through the unrolled iterations. Compute the paper-derived quantities exactly: A(x), b, c(x), and M. Estimate convergence by comparing against a very small-step reference trajectory over 32 independent random x-sequences, reporting RMS endpoint and path error. The first cheap experiment should use a small score model on CIFAR-10 or a synthetic stiff vector field containing a ReLU kink, comparing Euler, Heun, deterministic SDIRK, and randomized SDIRK at equal network evaluations. Success means lower trajectory error or better FID at fixed evaluations, together with fewer divergences in regimes where explicit Euler is unstable.