A library, written in Rust, for integrating ordinary differential equations. For now, this is relatively simple, but it does have key features that are needed for orbit propagation, ray tracing, and field line tracing:

Features

Explicit Runge-Kutta Methods (Non-Stiff Problems)

Method	Order	Stages	Dense Output	Best Use Case
BS3 (Bogacki-Shampine)	3(2)	4	3rd order	Moderate accuracy (rtol ~ 1e-4 to 1e-6)
DormandPrince45	5(4)	7	4th order	General purpose (rtol ~ 1e-6 to 1e-8)
Vern7 (Verner)	7(6)	10+6	7th order	High accuracy (rtol ~ 1e-8 to 1e-12)

Performance at 1e-10 tolerance:

Vern7: 2.7-8.8x faster than DP5
Vern7: 50x+ faster than BS3

See benchmark report for detailed performance analysis.

Other Features

Adaptive time stepping with PI controller
Callback events with zero-crossing detection
Dense output interpolation at any time point
Parameters in derivative and callback functions
Lazy computation of extra interpolation stages (Vern7)

Future Improvements

More algorithms
- Rosenbrock methods (for stiff problems)
- Tsit5
- Runge-Kutta Cash-Karp
Composite Algorithms
Automatic Stiffness Detection
Fixed Time Steps
Boolean callback eventing
Improved solution handling like DifferentialEquations.jl

To Use

For now, here is a simple example of using the propagator to solve a simple second-order system (the pendulum problem):

use nalgebra::Vector2;
use differential_equations::prelude::*;
use std::f64::consts::PI;

// Define the system (parameters, derivative, and initial state)
type Params = (f64, f64); // Gravity and Length of Pendulum
let params = (9.81, 1.0);

fn derivative(_t: f64, y: Vector2<f64>, p: &Params) -> Vector2<f64> {
    let &(g, l) = p;
    let theta = y[0];
    let d_theta = y[1];
    Vector2::new( d_theta, -(g/l) * theta.sin() )
}

let y0 = Vector2::new(0.0, PI/2.0);

// Set up the problem (ODE, Integrator, Controller, and Callbacks)
let ode = ODE::new(&derivative, 0.0, 6.3, y0, params);
let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-6);
let controller = PIController::default();

let value_too_high = Callback {
    event: &|t: f64, _y: Vector2<f64>, _p: &Params| { 5.0 - t },
    effect: &stop,
};

// Solve the problem
let mut problem = Problem::new(ode, dp45, controller).with_callback(value_too_high);
let solution = problem.solve();

// Can interpolate solutions to whatever you want
let _interpolated_answer = solution.interpolate(4.4);