Helpful context:

Most equations in algebra have numbers as solutions. The equation $x^2 - 4 = 0$ has the solution $x = 2$ (and $x = -2$). You find a number that makes the equation true.

But some equations have functions as solutions.

The equation $\frac{dy}{dt} = y$ asks something stranger: what function equals its own derivative? Whatever $y(t)$ is, its rate of change at every moment must equal its current value. If the function is currently at $3$, it must be growing at rate $3$. If it is at $100$, it must be growing at rate $100$.

Think about what kind of function would do this. If you double the current value, the derivative doubles too. The function grows in proportion to itself. That means it grows faster and faster as it gets larger - the more it has, the faster it grows. The answer is $y(t) = Ce^t$ for any constant $C$. Check: the derivative of $Ce^t$ is $Ce^t$. It equals itself exactly.

This equation, $\frac{dy}{dt} = y$, is an ordinary differential equation. And this type of equation - where a function’s rate of change depends on the function’s own value - describes an enormous range of real phenomena. Radioactive decay: the rate of decay is proportional to how much radioactive material remains. Compound interest: the rate of growth of a bank account is proportional to the current balance. Population growth (in simple models): the birth rate is proportional to the current population size. Cooling: the rate at which an object loses heat is proportional to the difference between its temperature and the ambient temperature (Newton’s law of cooling).

This post builds the theory of ordinary differential equations from the ground up, covering the key solution methods, the geometry of solutions, existence and uniqueness, and why ODEs appear throughout science and engineering.

Section 1: What Is a Differential Equation?

An ordinary differential equation (ODE) is an equation relating an unknown function $y(t)$ to one or more of its derivatives $y'(t), y''(t), \ldots$

The word “ordinary” distinguishes these from partial differential equations, which involve functions of multiple variables and their partial derivatives (temperature at each point in space and time, for instance). In an ODE, there is exactly one independent variable - here written $t$ for time, though it can represent any single variable.

The order of an ODE is the order of the highest derivative that appears:

$$\frac{dy}{dt} = ky \quad \text{(first order)}$$

$$\frac{d^2 y}{dt^2} + \omega^2 y = 0 \quad \text{(second order)}$$

$$y''' + y'' - y = e^t \quad \text{(third order)}$$

The solution is not a number. The solution is a function $y(t)$ that makes the equation true at every value of $t$.

Discomfort check. Why are solutions functions, not numbers?

Because the equation is a relationship that must hold at every moment. The equation $\frac{dy}{dt} = y$ is not asking “what number $y$ satisfies this?” It is asking “what function $y(t)$, when differentiated, gives back the same function at every time $t$?” The equation constrains $y$ across its entire domain simultaneously. A number has no derivative. Only a function does. The solution is an entire function.

Initial conditions. The general solution to $\frac{dy}{dt} = ky$ is $y(t) = Ce^{kt}$ for any constant $C$. The constant is undetermined until you specify the function’s value at some particular time. The equation $y(0) = y_0$ pins down $C = y_0$, giving the unique solution $y(t) = y_0 e^{kt}$. Such a specification is called an initial condition. An ODE together with initial conditions is an initial value problem (IVP).

For first-order ODEs, one initial condition determines the solution completely. For second-order ODEs, two initial conditions are needed - the value and the first derivative at some time $t_0$. In general, an $n$-th order ODE requires $n$ initial conditions.

Linear versus nonlinear. An ODE is linear if it is linear in $y$ and all its derivatives - no products like $y \cdot y'$ and no nonlinear terms like $\sin(y)$ or $y^2$. Linear ODEs have a rich algebraic theory. Nonlinear ODEs are much harder and often require qualitative or numerical methods.

Section 2: First-Order Separable ODEs

The simplest class of ODEs to solve directly is the separable first-order ODE:

$$\frac{dy}{dt} = f(y) \cdot g(t)$$

The right-hand side factors into a piece that depends only on $y$ and a piece that depends only on $t$. The strategy: put all $y$’s on one side and all $t$’s on the other, then integrate both sides.

$$\frac{dy}{f(y)} = g(t) dt \implies \int \frac{1}{f(y)} dy = \int g(t) dt$$

This algebraic separation is valid whenever $f(y) \neq 0$.

Example 1: Exponential growth and decay.

The ODE $\frac{dy}{dt} = ky$ where $k$ is a constant models growth ($k > 0$) or decay ($k < 0$). Separate:

$$\frac{dy}{y} = k dt$$

Integrate both sides:

$$\ln|y| = kt + C_0$$

Exponentiate: $|y| = e^{C_0} e^{kt}$. Writing $C = \pm e^{C_0}$ to absorb both the sign and the constant:

$$y(t) = Ce^{kt}$$

If $k > 0$: the function grows exponentially without bound. If $k < 0$: it decays toward zero. The rate is entirely controlled by $k$.

Application: half-life. The half-life $T_{1/2}$ is the time for $y$ to halve. Set $y(T_{1/2}) = y_0/2$:

$$y_0 e^{kT_{1/2}} = \frac{y_0}{2} \implies e^{kT_{1/2}} = \frac{1}{2} \implies T_{1/2} = \frac{-\ln 2}{k}$$

Carbon-14 has a half-life of approximately 5730 years, meaning $k = -\ln(2)/5730 \approx -0.000121$ per year. Radiocarbon dating measures the current ratio of C-14 to stable C-12, computes $y/y_0$, and solves $e^{kt} = y/y_0$ for $t$ - recovering the time since the organism died.

Example 2: Newton’s law of cooling.

A hot object cools in a room. Newton’s law says the rate of temperature change is proportional to the temperature difference between the object and the room:

$$\frac{dT}{dt} = -k(T - T_{\text{room}})$$

where $k > 0$ is a cooling constant. Let $u = T - T_{\text{room}}$ (the excess temperature). Then $\frac{du}{dt} = -ku$, which separates to $\frac{du}{u} = -k dt$, giving $u(t) = u_0 e^{-kt}$. Therefore:

$$T(t) = T_{\text{room}} + (T_0 - T_{\text{room}}) e^{-kt}$$

As $t \to \infty$, $T(t) \to T_{\text{room}}$: the object cools toward room temperature, approaching it asymptotically. This ODE is the mathematical model for the cooling curve of a computer chip, the warmth of a cup of coffee, and the temperature profile of a planet after a volcanic eruption.

Example 3: Logistic growth.

The pure exponential model $\frac{dy}{dt} = ky$ has a problem: it predicts unlimited growth. Real populations run into resource limits. The logistic equation introduces a carrying capacity $K$:

$$\frac{dy}{dt} = ky\left(1 - \frac{y}{K}\right)$$

When $y \ll K$: the factor $(1 - y/K) \approx 1$ and growth is approximately exponential. As $y \to K$: the factor approaches $0$ and growth slows. If $y$ somehow exceeds $K$: the factor is negative and the population decreases.

Separate:

$$\frac{dy}{y(1 - y/K)} = k dt$$

Partial fractions: $\frac{1}{y(1 - y/K)} = \frac{K}{y(K-y)} = \frac{1}{y} + \frac{1}{K - y}$. Integrate:

$$\ln|y| - \ln|K - y| = kt + C_0 \implies \frac{y}{K - y} = Ae^{kt}$$

Solve for $y$:

$$y(t) = \frac{K}{1 + Be^{-kt}}, \quad B = \frac{K - y_0}{y_0}$$

This is the logistic function: an S-shaped curve starting near $0$, growing rapidly through the middle, and asymptotically approaching $K$. It models epidemic curves (the “flatten the curve” shape), technology adoption, and is the direct ancestor of the neural network sigmoid $\sigma(x) = 1/(1 + e^{-x})$.

Where the sigmoid comes from. Set $K = 1$, $k = 1$, and let the initial time be $t = -\infty$ (so the curve has been growing forever). The formula gives $y(t) = 1/(1 + e^{-t}) = \sigma(t)$. The sigmoid activation function is the solution to the logistic ODE on an infinite time horizon. Its S-shape, asymptotes at $0$ and $1$, and inflection point at $t = 0$ are all consequences of the logistic dynamics.

Section 3: First-Order Linear ODEs

The first-order linear ODE in standard form:

$$\frac{dy}{dt} + p(t)y = q(t)$$

When $q(t) = 0$ this is separable. When $q(t) \neq 0$ we need the integrating factor method.

The idea. Find a function $\mu(t)$ such that multiplying both sides by $\mu(t)$ makes the left side a perfect derivative: $\frac{d}{dt}[\mu(t) y(t)]$.

By the product rule: $\frac{d}{dt}[\mu y] = \mu y' + \mu' y$. We want this to equal $\mu(y' + py) = \mu y' + \mu p y$. Matching requires $\mu' = \mu p$, which gives:

$$\mu(t) = e^{\int p(t) dt}$$

Multiply the ODE by $\mu(t)$:

$$\frac{d}{dt}[\mu(t) y(t)] = \mu(t) q(t)$$

Integrate both sides:

$$y(t) = \frac{1}{\mu(t)} \left(\int \mu(t) q(t) dt + C\right)$$

Worked example. Solve $\frac{dy}{dt} + \frac{y}{t} = t$ for $t > 0$ with $y(1) = 2$.

Here $p(t) = 1/t$ and $q(t) = t$. Integrating factor: $\mu = e^{\int (1/t) dt} = e^{\ln t} = t$.

Multiply through: $t y' + y = t^2$, so $\frac{d}{dt}[ty] = t^2$.

Integrate: $ty = \frac{t^3}{3} + C$, giving $y = \frac{t^2}{3} + \frac{C}{t}$.

Apply $y(1) = 2$: $2 = \frac{1}{3} + C$, so $C = \frac{5}{3}$.

The unique solution is $y(t) = \frac{t^2}{3} + \frac{5}{3t}$.

The homogeneous part $C/t$ decays as $t$ grows; the particular part $t^2/3$ grows. The initial condition determines which mixture you have.

Section 4: Second-Order Linear ODEs with Constant Coefficients

The equation:

$$ay'' + by' + cy = 0, \quad a \neq 0$$

This models spring-mass systems ($a$ = mass, $b$ = damping, $c$ = spring constant), RLC circuits, sound wave propagation, and a great deal else.

The characteristic equation. Guess $y = e^{rt}$. Then $y' = re^{rt}$ and $y'' = r^2 e^{rt}$. Substituting:

$$e^{rt}(ar^2 + br + c) = 0$$

Since $e^{rt} \neq 0$, the characteristic equation must hold:

$$ar^2 + br + c = 0$$

The discriminant $\Delta = b^2 - 4ac$ determines the nature of the roots, and hence the solution.

Case 1: Two distinct real roots ($\Delta > 0$).

Roots $r_1 \neq r_2$ via the quadratic formula. Two independent solutions $e^{r_1 t}$ and $e^{r_2 t}$. General solution:

$$y(t) = C_1 e^{r_1 t} + C_2 e^{r_2 t}$$

Discomfort check. Why do we need two arbitrary constants for a second-order ODE?

The solution space of a second-order linear ODE is a two-dimensional vector space. Just as $\mathbb{R}^2$ needs two basis vectors, this space needs two basis functions - two linearly independent solutions. Every solution is a linear combination of them. The two initial conditions ($y(t_0)$ and $y'(t_0)$) provide two equations that determine the two constants $C_1$ and $C_2$ uniquely. If you only had one constant, you could only satisfy one initial condition. Two constants, two equations, one unique solution.

Case 2: Repeated root ($\Delta = 0$).

One root $r = -b/(2a)$. The equation gives only one solution $e^{rt}$. A second independent solution is $te^{rt}$ (this can be verified by substitution, or derived by the technique of reduction of order). General solution:

$$y(t) = (C_1 + C_2 t)e^{rt}$$

Case 3: Complex roots ($\Delta < 0$).

Roots $r = \alpha \pm i\beta$ where $\alpha = -b/(2a)$ and $\beta = \sqrt{4ac - b^2}/(2a)$. Using Euler’s formula $e^{i\beta t} = \cos(\beta t) + i\sin(\beta t)$, the two complex solutions combine into two real solutions:

$$y(t) = e^{\alpha t}(C_1 \cos(\beta t) + C_2 \sin(\beta t))$$

Discomfort check. Complex roots produce real oscillatory solutions - where do the imaginary parts go?

The two roots $\alpha + i\beta$ and $\alpha - i\beta$ give complex-valued solutions $e^{(\alpha + i\beta)t}$ and $e^{(\alpha - i\beta)t}$. For a real ODE with real initial conditions, we want real solutions. Form the sum: $e^{(\alpha + i\beta)t} + e^{(\alpha - i\beta)t} = 2e^{\alpha t}\cos(\beta t)$. Form the difference divided by $2i$: $e^{\alpha t}\sin(\beta t)$. These are real, they are solutions (linear combinations of solutions are solutions for linear ODEs), and they are linearly independent. The imaginary parts cancel in forming these combinations. What remains: oscillation at frequency $\beta$ with amplitude modulated by $e^{\alpha t}$.

The spring-mass system in detail. The equation $mx'' + cx' + kx = 0$ with positive constants $m, c, k$ has characteristic roots:

$$r = \frac{-c \pm \sqrt{c^2 - 4mk}}{2m}$$

Note that $-c/(2m) < 0$ always, so $\alpha < 0$ always. The three regimes:

  • Underdamped ($c^2 < 4mk$): complex roots with $\alpha < 0$. Solution $x(t) = e^{\alpha t}(C_1 \cos\beta t + C_2 \sin\beta t)$ - oscillations with exponentially decaying amplitude. A plucked guitar string.

  • Overdamped ($c^2 > 4mk$): two distinct negative real roots. Solution decays exponentially without oscillation. A heavy door closer.

  • Critically damped ($c^2 = 4mk$): repeated negative root. Solution $(C_1 + C_2 t)e^{rt}$ - fastest return to equilibrium without oscillating. Engineered into car suspension and shock absorbers.

The transition between underdamped and overdamped at the critical damping condition $c^2 = 4mk$ is a fundamental design parameter in mechanical and electrical engineering.

Section 5: Systems of ODEs

Multiple coupled unknown functions:

$$\frac{d\mathbf{x}}{dt} = A\mathbf{x}$$

where $\mathbf{x}(t) \in \mathbb{R}^n$ is a vector and $A$ is an $n \times n$ constant matrix. This is the vector generalization of $dy/dt = ky$.

The matrix exponential. By analogy with the scalar solution $y(t) = e^{kt} y_0$, the solution is:

$$\mathbf{x}(t) = e^{At} \mathbf{x}(0)$$

where the matrix exponential is defined by the power series:

$$e^{At} = I + At + \frac{(At)^2}{2!} + \frac{(At)^3}{3!} + \cdots$$

This converges for every matrix $A$ and every $t \in \mathbb{R}$.

Computing via eigenvalues. If $A$ is diagonalizable with eigenvalues $\lambda_1, \ldots, \lambda_n$ and eigenvectors $v_1, \ldots, v_n$, write the initial condition as:

$$\mathbf{x}(0) = c_1 v_1 + c_2 v_2 + \cdots + c_n v_n$$

The solution decomposes into independent modes:

$$\mathbf{x}(t) = c_1 e^{\lambda_1 t} v_1 + c_2 e^{\lambda_2 t} v_2 + \cdots + c_n e^{\lambda_n t} v_n$$

Each eigenvalue determines the temporal behavior of its mode:

  • $\lambda_k < 0$ (negative real): mode decays exponentially to zero
  • $\lambda_k > 0$ (positive real): mode grows exponentially
  • $\lambda_k = \alpha \pm i\beta$ (complex pair): mode oscillates with amplitude $e^{\alpha t}$

The eigenvalues of $A$ determine the qualitative behavior of the entire system.

Stability. The system $\frac{d\mathbf{x}}{dt} = A\mathbf{x}$ is:

  • Asymptotically stable if all eigenvalues of $A$ have negative real part: every solution decays to $\mathbf{0}$.
  • Unstable if any eigenvalue has positive real part: some solutions grow without bound.
  • Neutrally stable if eigenvalues lie on the imaginary axis (no positive real parts): solutions oscillate without decaying.

This eigenvalue criterion is the foundation of linear stability analysis throughout engineering and applied mathematics.

Example. The predator-prey system linearized near an equilibrium. Two populations $x_1$ (prey) and $x_2$ (predators), with coupling matrix:

$$A = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}$$

Eigenvalues: characteristic polynomial $\lambda^2 + 1 = 0$, so $\lambda = \pm i$. Purely imaginary - neutrally stable. The solution oscillates periodically. This is the mathematical structure behind Lotka-Volterra cycles: predator and prey populations oscillate out of phase indefinitely, neither dying out nor exploding.

Section 6: Numerical Methods

Most ODEs arising in practice cannot be solved analytically. Numerical methods approximate the continuous solution using discrete steps.

Euler’s method. Given $\frac{dy}{dt} = f(t, y)$ with $y(t_0) = y_0$, choose step size $h$ and iterate:

$$y_{n+1} = y_n + h \cdot f(t_n, y_n), \quad t_{n+1} = t_n + h$$

The idea: at each point $(t_n, y_n)$, the ODE tells us the slope $f(t_n, y_n)$. Walk along that slope for a step of size $h$ to reach the next point.

Geometry. The ODE $\frac{dy}{dt} = f(t, y)$ defines a direction field (or slope field): at every point $(t, y)$ in the plane, draw a short line segment with slope $f(t, y)$. The solution curves are exactly the curves that are everywhere tangent to these segments. Euler’s method follows the direction field step by step - at each point, look at the slope, take a step in that direction, repeat.

Error. Over a fixed time interval $[t_0, T]$, Euler’s method has global error $O(h)$: halving the step size halves the error (approximately). This is modest accuracy - you need very small $h$ for precise results.

Runge-Kutta methods achieve higher accuracy by evaluating $f$ at multiple points within each step. The fourth-order Runge-Kutta (RK4) uses four evaluations per step:

$$k_1 = f(t_n, y_n)$$ $$k_2 = f\left(t_n + \tfrac{h}{2}, y_n + \tfrac{h}{2}k_1\right)$$ $$k_3 = f\left(t_n + \tfrac{h}{2}, y_n + \tfrac{h}{2}k_2\right)$$ $$k_4 = f(t_n + h, y_n + h k_3)$$ $$y_{n+1} = y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)$$

The global error of RK4 is $O(h^4)$: halving $h$ reduces the error by a factor of 16. This makes RK4 dramatically more efficient than Euler for the same accuracy.

Adaptive step-size control. Real ODE solvers (scipy’s solve_ivp with the DOPRI5 method, for instance) automatically adjust $h$ based on a local error estimate. Where the solution is smooth, they take large steps. Where it changes rapidly, they take small steps. The user specifies a tolerance, and the solver guarantees the error stays within it.

Stiff ODEs. An ODE is stiff when it has solution components evolving on vastly different timescales. The equation $y' = -1000y + 1$ has a fast transient (decaying with timescale $1/1000$) and a slow equilibrium at $y = 0.001$. Explicit methods like Euler and RK4 require step sizes $h < 2/1000 = 0.002$ for numerical stability - even when the interesting behavior is at timescale $1$. This is enormously wasteful. Implicit methods, which solve for $y_{n+1}$ using $f(t_{n+1}, y_{n+1})$ on the right-hand side, are stable for much larger steps. The tradeoff: each implicit step requires solving a nonlinear equation.

Section 7: Existence and Uniqueness

Before solving an ODE, it is worth knowing: does a solution exist? Is it unique?

Theorem (Picard-Lindelöf). Consider the initial value problem:

$$\frac{dy}{dt} = f(t, y), \quad y(t_0) = y_0$$

If $f$ is continuous in a rectangle $R = \{|t - t_0| \leq a, |y - y_0| \leq b\}$ and satisfies a Lipschitz condition in $y$ on $R$ - meaning there exists $L > 0$ such that:

$$|f(t, y_1) - f(t, y_2)| \leq L|y_1 - y_2|$$

for all $(t, y_1), (t, y_2) \in R$ - then there exists an interval $(t_0 - h, t_0 + h)$ on which the IVP has a unique solution.

The Lipschitz condition means $f$ cannot change too abruptly as $y$ varies. Most smooth functions satisfy it locally (if $f$ has a bounded partial derivative $\partial f / \partial y$ on $R$, it is Lipschitz with $L = \sup|\partial f / \partial y|$).

Proof sketch. Convert the ODE to an integral equation: $y(t) = y_0 + \int_{t_0}^t f(s, y(s)) ds$. Define the Picard iterates: $y_0(t) = y_0$, $y_{n+1}(t) = y_0 + \int_{t_0}^t f(s, y_n(s)) ds$. The Lipschitz condition makes the iteration a contraction on a function space (using the Banach fixed-point theorem), so the iterates converge to the unique fixed point - the solution.

Finite-time blowup. The theorem guarantees a local solution, but not a global one. Consider $\frac{dy}{dt} = y^2$, $y(0) = 1$. Separating: $\int y^{-2} dy = \int dt$ gives $-1/y = t + C$, so $y(t) = \frac{1}{1 - t}$. The solution blows up at $t = 1$ and cannot be extended beyond it. The Picard-Lindelöf theorem guarantees existence near $t = 0$ but makes no promise about $t = 1$.

Non-uniqueness without Lipschitz. The equation $\frac{dy}{dt} = y^{1/3}$, $y(0) = 0$ has the solution $y(t) = 0$ (check: $0 = 0$). But it also has $y(t) = (2t/3)^{3/2}$ for $t > 0$ (check: the derivative is $(2/3)^{3/2} \cdot (3/2) t^{1/2} = (2t/3)^{1/2}$, and $(y)^{1/3} = (2t/3)^{1/2}$). Two solutions from the same initial condition. The culprit: $f(y) = y^{1/3}$ has unbounded derivative near $y = 0$, violating the Lipschitz condition at the initial point.

Section 8: Phase Portraits

For autonomous systems $\frac{d\mathbf{x}}{dt} = F(\mathbf{x})$ (right-hand side does not explicitly depend on $t$), the qualitative behavior of all solutions can be visualized in the phase portrait: a plot of solution curves in the state space, with arrows showing the direction of motion.

Equilibria. Points where $F(\mathbf{x}^*) = 0$. At an equilibrium, $\frac{d\mathbf{x}}{dt} = 0$ - the system sits still. Perturbations either decay back (stable equilibrium) or grow away (unstable).

Stability via linearization. Near an equilibrium $\mathbf{x}^$, let $\mathbf{u} = \mathbf{x} - \mathbf{x}^$ be the small perturbation. Taylor-expand $F(\mathbf{x}) = F(\mathbf{x}^* + \mathbf{u}) \approx J\mathbf{u}$ where $J = \frac{\partial F}{\partial \mathbf{x}}|_{\mathbf{x}^*}$ is the Jacobian matrix. The linearized system $\frac{d\mathbf{u}}{dt} = J\mathbf{u}$ determines local stability: eigenvalues of $J$ with negative real parts attract, positive real parts repel.

Classification of 2D equilibria:

Eigenvalues of $J$ Type Stability
$\lambda_1, \lambda_2 < 0$, real Stable node Stable
$\lambda_1, \lambda_2 > 0$, real Unstable node Unstable
$\lambda_1 < 0 < \lambda_2$ Saddle point Unstable
$\alpha \pm i\beta$, $\alpha < 0$ Stable spiral Stable
$\alpha \pm i\beta$, $\alpha > 0$ Unstable spiral Unstable
$\pm i\beta$ Center Neutrally stable

The pendulum. The equation $\ddot{\theta} + (g/L)\sin\theta = 0$ is nonlinear. Written as a system with $\omega = \dot{\theta}$: the equilibrium at $(\theta, \omega) = (0, 0)$ (pendulum hanging down) linearizes to a center with purely imaginary eigenvalues - small oscillations. The equilibrium at $(\theta, \omega) = (\pi, 0)$ (pendulum balanced upright) linearizes to a saddle - unstable.

Limit cycles. Some nonlinear systems exhibit limit cycles - isolated closed orbits in the phase portrait. A trajectory that starts near the limit cycle spirals toward it as $t \to \infty$ (stable limit cycle) or away from it (unstable). Limit cycles are not possible in linear systems - they are a genuinely nonlinear phenomenon. The heart beats in a limit cycle: perturb the rhythm slightly and it returns to the regular oscillation.

The van der Pol oscillator $\ddot{x} - \mu(1 - x^2)\dot{x} + x = 0$ has a stable limit cycle for any $\mu > 0$. When $|x| < 1$, the damping coefficient $-\mu(1-x^2) < 0$ (negative damping - the system adds energy). When $|x| > 1$, the coefficient is positive (positive damping - the system removes energy). This self-regulating mechanism drives the trajectory to a specific amplitude, regardless of initial conditions.

Section 8b: A Complete Worked Example

Let us solve the initial value problem:

$$y'' - 5y' + 6y = e^{2t}, \quad y(0) = 1, \quad y'(0) = 0$$

from start to finish.

Step 1: Solve the homogeneous equation $y'' - 5y' + 6y = 0$.

Characteristic equation: $r^2 - 5r + 6 = 0$. Factor: $(r-2)(r-3) = 0$. Roots: $r_1 = 2$, $r_2 = 3$.

Homogeneous solution: $y_h = C_1 e^{2t} + C_2 e^{3t}$.

Step 2: Find a particular solution to $y'' - 5y' + 6y = e^{2t}$.

The natural guess $y_p = Ae^{2t}$ fails because $e^{2t}$ is already a homogeneous solution (it solves the homogeneous equation with the root $r = 2$). Multiply by $t$: guess $y_p = Ate^{2t}$.

Compute: $y_p' = Ae^{2t} + 2Ate^{2t} = A(1 + 2t)e^{2t}$. And $y_p'' = 2Ae^{2t} + 2Ae^{2t} + 4Ate^{2t} = A(4 + 4t)e^{2t}$.

Substitute into the ODE:

$$A(4 + 4t)e^{2t} - 5A(1 + 2t)e^{2t} + 6Ate^{2t} = e^{2t}$$

$$A e^{2t}[(4 + 4t) - 5(1 + 2t) + 6t] = e^{2t}$$

$$A e^{2t}[4 + 4t - 5 - 10t + 6t] = e^{2t}$$

$$A e^{2t}[-1 + 0 \cdot t] = e^{2t}$$

So $A(-1) = 1$, giving $A = -1$. Particular solution: $y_p = -te^{2t}$.

Step 3: General solution.

$$y(t) = y_h + y_p = C_1 e^{2t} + C_2 e^{3t} - te^{2t} = (C_1 - t)e^{2t} + C_2 e^{3t}$$

Step 4: Apply initial conditions.

$y(0) = C_1 + C_2 = 1$.

$y'(t) = (C_1 - t)2e^{2t} - e^{2t} + 3C_2 e^{3t} = (2C_1 - 2t - 1)e^{2t} + 3C_2 e^{3t}$.

$y'(0) = 2C_1 - 1 + 3C_2 = 0$.

System: $C_1 + C_2 = 1$ and $2C_1 + 3C_2 = 1$.

From the first: $C_1 = 1 - C_2$. Substitute: $2(1 - C_2) + 3C_2 = 1$, so $2 + C_2 = 1$, giving $C_2 = -1$ and $C_1 = 2$.

Final answer: $y(t) = (2 - t)e^{2t} - e^{3t}$.

Checking the answer: As $t \to \infty$, the term $-e^{3t}$ dominates (grows fastest), so $y(t) \to -\infty$. The system is unstable - both characteristic roots $r = 2$ and $r = 3$ are positive.

Section 9: The Structure of Linear Solutions

For linear homogeneous ODEs, the solution space has vector space structure.

Superposition. If $y_1(t)$ and $y_2(t)$ both satisfy $ay'' + by' + cy = 0$, then so does any linear combination $C_1 y_1 + C_2 y_2$. Proof: substitute $C_1 y_1 + C_2 y_2$ into the equation; by linearity of differentiation, you get $C_1(ay_1'' + by_1' + cy_1) + C_2(ay_2'' + by_2' + cy_2) = 0 + 0 = 0$.

The set of all solutions forms a vector space of dimension $n$ (for an $n$-th order linear ODE). Any set of $n$ linearly independent solutions is a fundamental set of solutions - a basis for this space.

The Wronskian. Two solutions $y_1, y_2$ are linearly independent if and only if their Wronskian is nonzero:

$$W(y_1, y_2)(t) = \det \begin{pmatrix} y_1(t) & y_2(t) \\ y_1'(t) & y_2'(t) \end{pmatrix} = y_1 y_2' - y_2 y_1' \neq 0$$

Abel’s theorem: if $y_1, y_2$ are solutions to $y'' + p(t)y' + q(t)y = 0$, then $W(t) = W(t_0) e^{-\int_{t_0}^t p(s) ds}$. The Wronskian is either identically zero (solutions are dependent) or never zero (solutions are independent).

Inhomogeneous equations. For $ay'' + by' + cy = g(t)$, the general solution is:

$$y(t) = y_h(t) + y_p(t)$$

where $y_h$ is the general homogeneous solution (with two free constants) and $y_p$ is any particular solution to the inhomogeneous equation.

Finding $y_p$: when $g(t)$ is a polynomial, exponential, sine, or cosine, the method of undetermined coefficients works - guess $y_p$ has the same form as $g$ and solve for coefficients. For general $g$, variation of parameters gives:

$$y_p(t) = -y_1(t) \int \frac{y_2 g}{W} dt + y_2(t) \int \frac{y_1 g}{W} dt$$

where $\{y_1, y_2\}$ is a fundamental set and $W$ is their Wronskian.

Section 10: Why ODEs Are Everywhere

The language of ODEs is the natural language of change. Any quantity whose rate of change depends on its current state satisfies an ODE.

Newton’s second law. $F = ma$ is $m\ddot{x} = F(x, \dot{x}, t)$. Every classical mechanics problem is an ODE - from projectile motion to celestial mechanics to molecular dynamics.

Electric circuits. The RLC circuit satisfies exactly the same equation as the spring-mass system: $L q'' + R q' + q/C = V(t)$, where $q(t)$ is charge, $L$ is inductance, $R$ is resistance, and $C$ is capacitance. The mathematical identity means every result proved for mechanical oscillators applies immediately to electrical oscillators.

Heat and diffusion. Newton’s law of cooling: $\frac{dT}{dt} = -k(T - T_{\text{ambient}})$. A separable ODE with solution $T(t) = T_{\text{ambient}} + (T_0 - T_{\text{ambient}})e^{-kt}$.

Chemical kinetics. The rate of a first-order reaction $A \to B$ satisfies $\frac{d[A]}{dt} = -k[A]$. Exponential decay of reactant concentration.

Neural dynamics. The Hodgkin-Huxley model of neuron action potentials is a system of four coupled nonlinear ODEs (membrane voltage plus three gating variables). Computational neuroscience is built on this dynamical systems framework.

Gradient flow. The continuous-time limit of gradient descent is the ODE $\frac{d\theta}{dt} = -\nabla L(\theta)$. Analysis of gradient descent - convergence rates, sensitivity to initialization - is analysis of this ODE. Neural ODEs take this connection seriously and make the neural network itself the right-hand side of an ODE.

Epidemic modeling. The SIR model for disease spread: $\frac{dS}{dt} = -\beta SI$, $\frac{dI}{dt} = \beta SI - \gamma I$, $\frac{dR}{dt} = \gamma I$, where $S$, $I$, $R$ are the sizes of susceptible, infected, and recovered populations. This system of three coupled ODEs has explained epidemic curves since Kermack and McKendrick (1927).

The deep point: once you know the rate of change - the right-hand side of the ODE - you know the entire future trajectory (given an initial condition and granted existence and uniqueness). ODEs convert local rules about rates into global predictions about futures.

Section 11: Reading a Solution

Given a solution $y(t)$, what does it tell you? Certain features are immediate.

Long-time behavior. The sign of the real parts of the characteristic roots determines whether the solution grows, decays, or oscillates. For $y'' + 3y' + 2y = 0$, the roots are $r = -1$ and $r = -2$ (both negative), so $y(t) = C_1 e^{-t} + C_2 e^{-2t} \to 0$ as $t \to \infty$.

Transients vs. steady state. For an inhomogeneous ODE with constant forcing, the homogeneous solution decays (if the system is stable) and the particular solution remains. The homogeneous part is the transient - the initial excitement that fades. The particular part is the steady state - the long-run behavior.

Resonance. If the forcing frequency equals a natural frequency of the system, the particular solution grows. For $y'' + \omega^2 y = \cos(\omega t)$ (forcing at exactly the natural frequency), the undetermined coefficients method must be modified: the usual guess $A\cos(\omega t) + B\sin(\omega t)$ is already a homogeneous solution. The correct particular solution is $y_p = \frac{t}{2\omega}\sin(\omega t)$ - amplitude grows linearly in time. This is resonance, and it is why soldiers break step when crossing bridges.

Energy in oscillating systems. For the simple harmonic oscillator $\ddot{x} + \omega^2 x = 0$, the total mechanical energy is $E = \frac{1}{2}\dot{x}^2 + \frac{1}{2}\omega^2 x^2$. Differentiating and using the ODE: $\dot{E} = \dot{x}\ddot{x} + \omega^2 x \dot{x} = \dot{x}(\ddot{x} + \omega^2 x) = 0$. The energy is conserved. The solution curves in the phase plane $(x, \dot{x})$ are level sets of $E$ - ellipses centered at the origin. This conserved quantity is a first integral of the system, and it constrains where solutions can go.

Long-time behavior and stability redux. Whether a solution stays bounded, grows, or decays tells you about the stability of the corresponding equilibrium. For the system $y'' + 3y' + 2y = 0$, the roots are $r = -1, -2$. Both negative: every solution decays to zero. For $y'' - y = 0$, the roots are $r = \pm 1$. One root is positive: solutions starting away from the origin grow exponentially. The sign of the rightmost characteristic root is the key indicator.

When all characteristic roots have negative real parts, the ODE is called exponentially stable: every solution decays like $Ce^{-\mu t}$ for some $\mu > 0$. This is the property you want in control systems (perturbations die out) and trained neural networks (the training gradient ODE stays bounded).

Summary

Concept Key idea
ODE Equation relating a function to its derivatives; solution is a function
Order Highest derivative appearing; determines dimension of solution space
Separable ODE $dy/f(y) = g(t) dt$; integrate both sides
Integrating factor $\mu = e^{\int p dt}$; converts $y' + py = q$ to $d[\mu y]/dt = \mu q$
Characteristic equation $ar^2 + br + c = 0$; roots determine homogeneous solution form
Distinct real roots $y = C_1 e^{r_1 t} + C_2 e^{r_2 t}$
Repeated root $y = (C_1 + C_2 t)e^{rt}$
Complex roots $\alpha \pm i\beta$ $y = e^{\alpha t}(C_1 \cos\beta t + C_2 \sin\beta t)$
System $\dot{\mathbf{x}} = A\mathbf{x}$ Decomposes into eigenvalue modes $c_k e^{\lambda_k t} v_k$
Stability All eigenvalues with $\text{Re}(\lambda) < 0$ means stable
Euler method $y_{n+1} = y_n + hf(t_n, y_n)$; global error $O(h)$
RK4 Four evaluations per step; global error $O(h^4)$
Picard-Lindelöf Lipschitz condition guarantees local existence and uniqueness
Superposition For linear ODEs, linear combinations of solutions are solutions
Wronskian $W = y_1 y_2' - y_2 y_1'$; nonzero iff solutions are independent

Read Next: