Solving Differential Equations Using MathCore Rust

MathCore’s differential module provides a complete toolkit for solving ordinary and partial differential equations numerically. ODEs are expressed as symbolic Expr values — you write the right-hand side as a MathCore expression string, and the solver evaluates it step by step. PDE solvers use finite-difference schemes and operate on closures for boundary and initial conditions. All solvers return structured result types that make downstream plotting and analysis straightforward.

Import

use mathcore::differential::{DifferentialEquations, PDESolver};
use mathcore::MathCore;

Result Types

`ODESolution`

All ODE solvers return an ODESolution struct:

#[derive(Debug, Clone)]
pub struct ODESolution {
    pub t: Vec<f64>,       // time (or independent variable) values
    pub y: Vec<Vec<f64>>,  // solution values; y[i] is the state vector at t[i]
}

For scalar (first-order, single equation) solvers, y[i] is a one-element vector [y_i]. For systems of ODEs, y[i] contains one entry per equation.

`BoundaryCondition`

#[derive(Debug, Clone)]
pub enum BoundaryCondition {
    InitialValue {
        t0: f64,
        y0: Vec<f64>,
    },
    BoundaryValue {
        ta: f64,
        tb: f64,
        ya: Vec<f64>,
        yb: Vec<f64>,
    },
}

ODE Solvers

First-Order ODE — RK4

`DifferentialEquations::solve_ode_first_order(f, t_var, y_var, initial_condition, t_final, steps) -> Result<ODESolution, MathError>`

Solves dy/dt = f(t, y) using the 4th-order Runge-Kutta method. This is the recommended solver for smooth, non-stiff problems.

Parameter	Type	Description
`f`	`&Expr`	Right-hand side expression; may reference `t_var` and `y_var`
`t_var`	`&str`	Name of the independent variable (e.g. `"t"`)
`y_var`	`&str`	Name of the dependent variable (e.g. `"y"`)
`initial_condition`	`(f64, f64)`	`(t₀, y₀)` — initial point
`t_final`	`f64`	End of the integration interval
`steps`	`usize`	Number of uniform steps; step size h = (t_final − t₀) / steps

The RK4 update formula applied at each step:

k1 = f(t,       y)
k2 = f(t + h/2, y + h·k1/2)
k3 = f(t + h/2, y + h·k2/2)
k4 = f(t + h,   y + h·k3)
y_new = y + h·(k1 + 2·k2 + 2·k3 + k4) / 6

Exponential decay — dy/dt = −y, y(0) = 1 (analytical: y = exp(-t)):

use mathcore::differential::DifferentialEquations;
use mathcore::MathCore;

let f = MathCore::parse("-y").unwrap();

let solution = DifferentialEquations::solve_ode_first_order(
    &f,
    "t",           // independent variable
    "y",           // dependent variable
    (0.0, 1.0),   // initial: t=0, y=1
    5.0,           // solve to t=5
    100,           // 100 steps
).unwrap();

println!("t values: {:?}", &solution.t[..5]);
println!("y values: {:?}", &solution.y[..5]);

// Verify: y(1) ≈ e^(-1) ≈ 0.3679
let final_y = solution.y.last().unwrap()[0];
println!("y(5) ≈ {:.6}", final_y); // ≈ 0.006738

Logistic growth — real-world example from examples/scientific_computing.rs:

use mathcore::differential::DifferentialEquations;
use mathcore::parser::Parser;

// dP/dt = r·P·(1 - P/K),  r=0.5, K=1000
let expr = Parser::parse("0.5*P*(1 - P/1000)").unwrap();

let solution = DifferentialEquations::solve_ode_first_order(
    &expr,
    "t",
    "P",
    (0.0, 10.0), // initial population 10
    20.0,         // 20 time units
    200,
).unwrap();

for i in (0..solution.t.len()).step_by(20) {
    println!("t={:.1}  P={:.1}", solution.t[i], solution.y[i][0]);
}

First-Order ODE — Euler’s Method

`DifferentialEquations::euler_method(f, t_var, y_var, initial_condition, t_final, steps) -> Result<ODESolution, MathError>`

Simple forward Euler integration: y_ = y_n + h·f(t_n, y_n). Same signature as solve_ode_first_order.

use mathcore::differential::DifferentialEquations;
use mathcore::MathCore;

let f = MathCore::parse("-y").unwrap();

let solution = DifferentialEquations::euler_method(
    &f, "t", "y",
    (0.0, 1.0),
    5.0,
    100,
).unwrap();

Euler’s method has first-order accuracy (local error O(h²), global error O(h)). RK4 achieves fourth-order accuracy with the same step count. Prefer solve_ode_first_order unless you have a specific reason to use Euler — such as debugging or demonstrating the accuracy difference.

Stiff ODE Solver

`DifferentialEquations::solve_stiff_ode(f, jacobian, t_var, y_var, initial_condition, t_final, steps, tolerance) -> Result<ODESolution, MathError>`

Solves stiff problems using the implicit Backward Euler method with Newton iteration at each step. Stiff ODEs have rapidly decaying components that force explicit methods to use very small step sizes.

Parameter	Type	Description
`jacobian`	`Option<&Expr>`	Analytic ∂f/∂y; if `None`, a numerical Jacobian (finite differences) is used
`tolerance`	`f64`	Newton iteration convergence threshold

use mathcore::differential::DifferentialEquations;
use mathcore::MathCore;

// Stiff equation: dy/dt = -1000*(y - cos(t))
let f = MathCore::parse("-1000*(y - cos(t))").unwrap();

let solution = DifferentialEquations::solve_stiff_ode(
    &f,
    None,          // use numerical Jacobian
    "t", "y",
    (0.0, 1.0),
    2.0,
    200,
    1e-8,          // Newton tolerance
).unwrap();

Providing an analytic Jacobian expression speeds up the Newton iteration significantly for stiff problems, because it avoids the extra function evaluation needed for the numerical approximation.

Second-Order ODE

`DifferentialEquations::solve_ode_second_order(p, q, r, x_var, initial_conditions, x_final, steps) -> Result<ODESolution, MathError>`

Solves the linear second-order equation y″ + p(x)·y′ + q(x)·y = r(x) by converting it to a first-order system internally. The returned ODESolution stores both y and y′ at each step: solution.y[i] = [y_i, y'_i].

Parameter	Type	Description
`p`, `q`, `r`	`&Expr`	Coefficient functions of the independent variable
`initial_conditions`	`(f64, f64, f64)`	`(x₀, y(x₀), y′(x₀))`

use mathcore::differential::DifferentialEquations;
use mathcore::MathCore;

// Simple harmonic oscillator: y'' + y = 0
// p=0, q=1, r=0; y(0)=1, y'(0)=0  =>  y = cos(x)
let p = MathCore::parse("0").unwrap();
let q = MathCore::parse("1").unwrap();
let r = MathCore::parse("0").unwrap();

let solution = DifferentialEquations::solve_ode_second_order(
    &p, &q, &r,
    "x",
    (0.0, 1.0, 0.0), // x0=0, y0=1, y'0=0
    6.28318,           // one full period (2π)
    1000,
).unwrap();

// solution.y[i][0] = y(x_i), solution.y[i][1] = y'(x_i)
println!("y(π) ≈ {:.4}", solution.y[500][0]); // ≈ -1

System of First-Order ODEs

`DifferentialEquations::solve_ode_system(functions, t_var, y_vars, initial_conditions, t_final, steps) -> Result<ODESolution, MathError>`

Solves dy⃗/dt = F⃗(t, y⃗) using RK4 applied to all equations simultaneously. Each equation in functions may reference t_var and any variable name in y_vars. Nonlinear pendulum — from examples/scientific_computing.rs:

use mathcore::differential::DifferentialEquations;
use mathcore::parser::Parser;

// θ' = ω,  ω' = -(g/L)·sin(θ)
let theta_dot = Parser::parse("omega").unwrap();
let omega_dot = Parser::parse("-9.81*sin(theta)").unwrap();

let solution = DifferentialEquations::solve_ode_system(
    &[theta_dot, omega_dot],
    "t",
    &["theta".to_string(), "omega".to_string()],
    (0.0, vec![1.0, 0.0]), // θ(0)=1 rad, ω(0)=0
    10.0,                   // 10 seconds
    1000,
).unwrap();

println!("Time\tAngle\t\tAngular Velocity");
for i in (0..solution.t.len()).step_by(100) {
    println!(
        "{:.1}\t{:.4}\t\t{:.4}",
        solution.t[i],
        solution.y[i][0], // theta
        solution.y[i][1], // omega
    );
}

Analytical Solution for Linear ODEs with Constant Coefficients

`DifferentialEquations::solve_linear_constant_coeff(coefficients: &[f64], initial_conditions: &[f64]) -> Result<Expr, MathError>`

Returns a symbolic general solution for an nth-order linear ODE with constant coefficients. coefficients must be provided from highest to lowest degree: [a_n, a_{n-1}, …, a_1, a_0]. Real roots produce exponential terms Cᵢ·exp(rᵢ·t) and complex conjugate pairs produce exp(α·t)·(C₁·cos(β·t) + C₂·sin(β·t)) terms.

use mathcore::differential::DifferentialEquations;

// y'' - 3y' + 2y = 0  =>  coefficients [1, -3, 2]
// Roots: r=1, r=2  =>  solution: C1*e^t + C2*e^(2t)
let solution = DifferentialEquations::solve_linear_constant_coeff(
    &[1.0, -3.0, 2.0],  // coefficients
    &[1.0, 0.0],         // initial conditions (y(0), y'(0))
).unwrap();

println!("General solution: {}", solution);

PDE Solvers

PDESolver uses finite-difference methods on uniform grids.

Heat Equation

`PDESolver::solve_heat_equation(alpha, initial_condition, boundary_left, boundary_right, x_range, t_final, nx, nt) -> Result<Vec<Vec<f64>>, MathError>`

Solves ∂u/∂t = α·∂²u/∂x² using the explicit finite-difference scheme. Returns a 2D grid u[t_index][x_index].

The explicit scheme is only stable when α·Δt/Δx² ≤ 0.5 (Von Neumann stability condition). If this ratio exceeds 0.5, the solver returns Err(MathError::InvalidOperation). Reduce t_final / nt or increase nx to satisfy the condition.

use mathcore::differential::PDESolver;

let initial_temp = |x: f64| {
    if (x - 0.5).abs() < 0.1 { 100.0 } else { 20.0 }
};

let grid = PDESolver::solve_heat_equation(
    0.01,           // thermal diffusivity α
    &initial_temp,
    20.0,           // left boundary temperature
    20.0,           // right boundary temperature
    (0.0, 1.0),     // x ∈ [0, 1]
    0.5,            // solve to t=0.5
    50,             // 50 spatial points
    100,            // 100 time steps
).unwrap();

println!("Initial max temp: {:.1}", grid[0].iter().cloned().fold(f64::NEG_INFINITY, f64::max));
println!("Final max temp:   {:.1}", grid[99].iter().cloned().fold(f64::NEG_INFINITY, f64::max));

Wave Equation

`PDESolver::solve_wave_equation(c, initial_position, initial_velocity, x_range, t_final, nx, nt) -> Result<Vec<Vec<f64>>, MathError>`

Solves ∂²u/∂t² = c²·∂²u/∂x² using the explicit Courant–Friedrichs–Lewy (CFL) scheme.

The CFL condition requires c·Δt/Δx ≤ 1. The solver returns an error if this is violated. Increase nt (more time steps) or decrease nx (coarser spatial grid) to satisfy it.

use mathcore::differential::PDESolver;

let initial_position = |x: f64| {
    let (center, width) = (0.5, 0.05);
    (-(((x - center) / width).powi(2))).exp()
};
let initial_velocity = |_: f64| 0.0;

let grid = PDESolver::solve_wave_equation(
    0.5,                  // wave speed c
    &initial_position,
    &initial_velocity,
    (0.0, 1.0),           // x ∈ [0, 1]
    1.0,                  // t ∈ [0, 1]
    100,                  // spatial points
    200,                  // time steps
).unwrap();

Laplace Equation

`PDESolver::solve_laplace_equation(boundary_conditions, x_range, y_range, nx, ny, tolerance, max_iterations) -> Result<Vec<Vec<f64>>, MathError>`

Solves ∇²u = 0 (steady-state heat distribution, electrostatics, etc.) using Gauss-Seidel iteration on a 2D grid. boundary_conditions is a closure (x, y) -> Option<f64>; return Some(value) on boundary points and None for interior points.

use mathcore::differential::PDESolver;

let bcs = |x: f64, y: f64| -> Option<f64> {
    if x == 0.0 || x == 1.0 || y == 0.0 { Some(0.0) }
    else if y == 1.0 { Some(100.0) }  // hot top wall
    else { None }
};

let grid = PDESolver::solve_laplace_equation(
    &bcs,
    (0.0, 1.0),   // x range
    (0.0, 1.0),   // y range
    50,            // nx
    50,            // ny
    1e-6,          // convergence tolerance
    10_000,        // max iterations
).unwrap();

println!("Center temperature: {:.2}", grid[25][25]);

RK4 vs Euler Accuracy

The choice of solver matters. For the same step size h, RK4 makes four function evaluations per step and achieves global error O(h⁴), while Euler makes one evaluation and produces O(h) error.

Solver	Evaluations/step	Global error	Best for
`euler_method`	1	O(h)	Simple demos, very cheap f
`solve_ode_first_order` (RK4)	4	O(h⁴)	Most smooth problems
`solve_stiff_ode` (Backward Euler)	~10 (Newton)	O(h)	Stiff problems

For a given accuracy target, RK4 typically allows a step size ~10× larger than Euler, meaning far fewer steps for the same precision. Unless your right-hand side f is extremely expensive to evaluate, RK4 is the better choice.

Get Started

Core Concepts

Calculus

Algebra & Solving

Advanced Features

Guides

Solving Differential Equations Using MathCore Rust

Import

Result Types

`ODESolution`

`BoundaryCondition`

ODE Solvers

First-Order ODE — RK4

`DifferentialEquations::solve_ode_first_order(f, t_var, y_var, initial_condition, t_final, steps) -> Result<ODESolution, MathError>`

First-Order ODE — Euler’s Method

`DifferentialEquations::euler_method(f, t_var, y_var, initial_condition, t_final, steps) -> Result<ODESolution, MathError>`

Stiff ODE Solver

`DifferentialEquations::solve_stiff_ode(f, jacobian, t_var, y_var, initial_condition, t_final, steps, tolerance) -> Result<ODESolution, MathError>`

Second-Order ODE

`DifferentialEquations::solve_ode_second_order(p, q, r, x_var, initial_conditions, x_final, steps) -> Result<ODESolution, MathError>`

System of First-Order ODEs

`DifferentialEquations::solve_ode_system(functions, t_var, y_vars, initial_conditions, t_final, steps) -> Result<ODESolution, MathError>`

Analytical Solution for Linear ODEs with Constant Coefficients

`DifferentialEquations::solve_linear_constant_coeff(coefficients: &[f64], initial_conditions: &[f64]) -> Result<Expr, MathError>`

PDE Solvers

Heat Equation

`PDESolver::solve_heat_equation(alpha, initial_condition, boundary_left, boundary_right, x_range, t_final, nx, nt) -> Result<Vec<Vec<f64>>, MathError>`

Wave Equation

`PDESolver::solve_wave_equation(c, initial_position, initial_velocity, x_range, t_final, nx, nt) -> Result<Vec<Vec<f64>>, MathError>`

Laplace Equation

`PDESolver::solve_laplace_equation(boundary_conditions, x_range, y_range, nx, ny, tolerance, max_iterations) -> Result<Vec<Vec<f64>>, MathError>`

RK4 vs Euler Accuracy

Build docs developers (and LLMs) love

Get Started

Core Concepts

Calculus

Algebra & Solving

Advanced Features

Guides

Documentation Index

​Import

​Result Types

​ODESolution

​BoundaryCondition

​ODE Solvers

​First-Order ODE — RK4

​DifferentialEquations::solve_ode_first_order(f, t_var, y_var, initial_condition, t_final, steps) -> Result<ODESolution, MathError>

​First-Order ODE — Euler’s Method

​DifferentialEquations::euler_method(f, t_var, y_var, initial_condition, t_final, steps) -> Result<ODESolution, MathError>

​Stiff ODE Solver

​DifferentialEquations::solve_stiff_ode(f, jacobian, t_var, y_var, initial_condition, t_final, steps, tolerance) -> Result<ODESolution, MathError>

​Second-Order ODE

​DifferentialEquations::solve_ode_second_order(p, q, r, x_var, initial_conditions, x_final, steps) -> Result<ODESolution, MathError>

​System of First-Order ODEs

​DifferentialEquations::solve_ode_system(functions, t_var, y_vars, initial_conditions, t_final, steps) -> Result<ODESolution, MathError>

​Analytical Solution for Linear ODEs with Constant Coefficients

​DifferentialEquations::solve_linear_constant_coeff(coefficients: &[f64], initial_conditions: &[f64]) -> Result<Expr, MathError>

​PDE Solvers

​Heat Equation

​PDESolver::solve_heat_equation(alpha, initial_condition, boundary_left, boundary_right, x_range, t_final, nx, nt) -> Result<Vec<Vec<f64>>, MathError>

​Wave Equation

​PDESolver::solve_wave_equation(c, initial_position, initial_velocity, x_range, t_final, nx, nt) -> Result<Vec<Vec<f64>>, MathError>

​Laplace Equation

​PDESolver::solve_laplace_equation(boundary_conditions, x_range, y_range, nx, ny, tolerance, max_iterations) -> Result<Vec<Vec<f64>>, MathError>

​RK4 vs Euler Accuracy

Build docs developers (and LLMs) love

Import

Result Types

`ODESolution`

`BoundaryCondition`

ODE Solvers

First-Order ODE — RK4

`DifferentialEquations::solve_ode_first_order(f, t_var, y_var, initial_condition, t_final, steps) -> Result<ODESolution, MathError>`

First-Order ODE — Euler’s Method

`DifferentialEquations::euler_method(f, t_var, y_var, initial_condition, t_final, steps) -> Result<ODESolution, MathError>`

Stiff ODE Solver

`DifferentialEquations::solve_stiff_ode(f, jacobian, t_var, y_var, initial_condition, t_final, steps, tolerance) -> Result<ODESolution, MathError>`

Second-Order ODE

`DifferentialEquations::solve_ode_second_order(p, q, r, x_var, initial_conditions, x_final, steps) -> Result<ODESolution, MathError>`

System of First-Order ODEs

`DifferentialEquations::solve_ode_system(functions, t_var, y_vars, initial_conditions, t_final, steps) -> Result<ODESolution, MathError>`

Analytical Solution for Linear ODEs with Constant Coefficients

`DifferentialEquations::solve_linear_constant_coeff(coefficients: &[f64], initial_conditions: &[f64]) -> Result<Expr, MathError>`

PDE Solvers

Heat Equation

`PDESolver::solve_heat_equation(alpha, initial_condition, boundary_left, boundary_right, x_range, t_final, nx, nt) -> Result<Vec<Vec<f64>>, MathError>`

Wave Equation

`PDESolver::solve_wave_equation(c, initial_position, initial_velocity, x_range, t_final, nx, nt) -> Result<Vec<Vec<f64>>, MathError>`

Laplace Equation

`PDESolver::solve_laplace_equation(boundary_conditions, x_range, y_range, nx, ny, tolerance, max_iterations) -> Result<Vec<Vec<f64>>, MathError>`

RK4 vs Euler Accuracy