Differential Equation Solver

A first-order linear ODE with constant coefficients takes the form dy/dx + py = q, where p is a real constant multiplying the unknown function y(x), and q is a constant forcing term (or source). The goal is to find a function y(x) that satisfies the equation for all x, along with a specified initial condition y(x₀) = y₀. These problems arise naturally whenever the rate of change of a quantity is proportional to its current value plus an external input — a remarkably common pattern in the physical and social sciences.

The integrating factor method is the standard approach for solving these equations analytically. The integrating factor is μ(x) = e^(px). Multiplying both sides of the ODE by μ(x) converts the left side into the exact derivative d/dx[μ(x)·y(x)], making the equation directly integrable. For constant coefficients, this yields the general solution y(x) = q/p + C·e^(−px), where the equilibrium value q/p is the particular solution and C·e^(−px) is the complementary (homogeneous) solution. When p = 0, the equation reduces to dy/dx = q with the simple linear solution y(x) = y₀ + q(x − x₀). The integration constant C is determined uniquely by substituting the initial condition y(x₀) = y₀, giving C = (y₀ − q/p)·e^(px₀).

Practical applications of this solver span many disciplines. In electrical engineering, an RL circuit satisfies dI/dt + (R/L)I = V/L — exactly this form — and the solution describes how current builds up toward a steady state. In thermodynamics, Newton's law of cooling gives dT/dt + kT = kT_env. In pharmacokinetics, drug concentration in the bloodstream follows dC/dt + kC = dose_rate. In finance, continuously compounded savings accounts with constant deposits obey the same structure. Understanding the equilibrium value q/p is especially useful: it represents the long-run behaviour of the system as x → ∞ (when p > 0), and the term C·e^(−px) describes the transient decay toward that equilibrium.

This solver handles only first-order linear ODEs with constant coefficients p and q. It does not support variable-coefficient equations such as dy/dx + P(x)y = Q(x) with non-constant P or Q, nor does it handle nonlinear ODEs (e.g., Bernoulli or Riccati equations), systems of ODEs, or higher-order equations. When p = 0, the formula reduces to a simple linear ramp and the equilibrium output is undefined — this case is handled separately in the calculation. Very large values of |p·x| may cause numerical overflow or underflow in the exponential terms; for such cases, consider rescaling variables. The solver assumes x is a real variable and does not perform complex arithmetic. Initial conditions with very large magnitudes combined with large |p| values can produce results that are numerically indistinguishable from the equilibrium value due to floating-point precision limits. For stiff equations or non-constant coefficients, numerical methods such as Euler's method, Runge-Kutta (RK4), or dedicated ODE software (MATLAB ode45, Python scipy.integrate.solve_ivp) should be used instead.