feat: Update BDF integrator to match SciPy implementation

- Replace history-based approach with SciPy's differences array (D)
- Implement SciPy BDF coefficients (gamma, alpha, error_const, kappa)
- Add compute_R and change_D functions for step size changes
- Implement proper Newton system solver with LU decomposition approach
- Update step size control with MIN_FACTOR, MAX_FACTOR, safety factors
- Implement sophisticated order selection strategy based on error norms
- Update error estimation using SciPy's error constants
- All tests pass in <3ms, well under 1 second target
- Maintains numerical accuracy while improving performance

Author: WenyinWei