Benchmarked custom PDE solver against analytical solutions: χ² = 0.005. What else should I test?

I've developed a PDE solver in Mathematica using hierarchical subdivision and want to ensure it's actually working correctly before applying to harder problems. Current benchmarks: Hydrogen atom ground state: E = -13.5989 eV (exact: -13.6000 eV, error 0.008%) 2D harmonic oscillator: energies match within 0.01% Wave packet dispersion: agrees with analytical Gaussian spreading Convergence study: Three grid resolutions show consistent error reduction. RMSE ratio (coarse/fine) ≈ 0.98-0.99, suggesting near-asymptotic regime. Max residual: |iℏ∂ψ/∂t - Ĥψ| ~ 3.08-3.09 across all resolutions (should decrease with refinement - why isn't it?) Question: What additional validation tests should I run? I want to be confident this isn't producing correct answers for wrong reasons. Are there standard test suites for quantum PDE solvers? Available: Can share notebooks if anyone wants to spot numerical issues I'm missing. Particularly interested in feedback on whether my time-stepping is actually unitary or just approximately so.