From a stack of equations to one vector equation
So far every equation in this rung has had one unknown function. But the world rarely offers them one at a time. Two tanks exchanging brine, two masses joined by springs, a predator and its prey, the currents in a two-loop circuit — in each, several quantities change together, and the rate of each depends on all the others. That coupling is the whole point: you cannot solve one without the rest. A system of linear ODEs packages this honestly. Write the unknowns as a vector x(t) = [x_1(t); x_2(t); ...; x_n(t)], and the rule of motion becomes dx/dt = A x, where A is a constant n-by-n matrix of coefficients. We will stay with constant A, the case that is both solvable and everywhere.
Two facts make this more than tidy notation. First, a single higher-order equation is secretly a system in disguise. Take the damped oscillator u'' + b u' + k u = 0. Name the velocity v = u', so that u' = v and v' = -k u - b v; in vector form, with x = [u; v], you get dx/dt = A x where A = [0, 1; -k, -b]. Every n-th order linear equation flattens this way into a first-order system of size n by introducing the derivatives as new variables. Systems are not a separate, harder topic — they are the natural language in which all of linear ODE theory finally speaks.
Second, the whole structural theory from the first guide in this rung survives intact, just dressed in vectors. The system dx/dt = A x is homogeneous and linear, so the superposition principle still holds: add two solution vectors, scale them, and you stay inside the solution set. The solutions of an n-dimensional system form an n-dimensional vector space, and we need n independent solution vectors to span it. Picture x(t) not as numbers but as a moving point — a particle traced out in the n-dimensional state space, with A x the vector field of velocities telling the point where to head at each location. Solving the system means finding how that point flows.
The eigenvalue method: guess a straight-line solution
How do we crack dx/dt = A x? Steal the idea that solved the single equation. For one scalar equation x' = a x, the answer was x = e^{at}: an exponential, because differentiating it just multiplies by a. Try the vector analogue — guess that the system has a solution which keeps the same direction forever and only grows or shrinks in time: x(t) = e^{lambda t} v, where v is a fixed nonzero vector and lambda a number. Differentiate: dx/dt = lambda e^{lambda t} v. Plug into the equation: lambda e^{lambda t} v = A (e^{lambda t} v) = e^{lambda t} A v. Cancel the never-zero scalar e^{lambda t}, and the calculus has completely vanished, leaving a pure algebra problem: A v = lambda v.
That equation A v = lambda v is the eigenvalue problem of linear algebra. A vector v whose direction A leaves unchanged — A merely stretches it by the factor lambda — is an eigenvector, and lambda is its eigenvalue. This is the eigenvalue method: each eigenpair (lambda, v) of A hands you one straight-line solution e^{lambda t} v, a trajectory that rides forever along the line through v, expanding if lambda > 0, collapsing toward the origin if lambda < 0. The eigenvectors are the special directions along which the coupled tangle of equations decouples into simple, independent exponential growth — the system's natural axes.
To find the eigenvalues, rewrite A v = lambda v as (A - lambda I) v = 0, where I is the identity matrix. A nonzero vector v can satisfy this only if the matrix A - lambda I crushes it to zero, which happens exactly when that matrix is singular — when det(A - lambda I) = 0. This determinant is a polynomial in lambda, the characteristic polynomial, and its roots are the eigenvalues. Notice the lovely echo: for the single damped oscillator you solved a characteristic equation for the exponents r; here the very same kind of polynomial, born now from a determinant, hands you the exponents lambda. The single equation was always just the 1-by-1 (or companion-matrix) shadow of this.
Solve dx/dt = A x with A = [ 1, 1 ; 4, 1 ]
Characteristic polynomial:
det(A - lambda I) = | 1-lambda 1 |
| 4 1-lambda |
= (1-lambda)^2 - 4 = lambda^2 - 2 lambda - 3 = (lambda-3)(lambda+1)
=> lambda = 3, lambda = -1
Eigenvectors:
lambda = 3 : (A-3I)v=0 -> -2 v1 + v2 = 0 -> v = [1; 2]
lambda = -1: (A+I)v=0 -> 2 v1 + v2 = 0 -> v = [1; -2]
General solution (superpose the two straight-line solutions):
x(t) = c1 e^{3t} [1; 2] + c2 e^{-1 t} [1; -2]The fundamental matrix: bundling the solutions
For an n-dimensional system you collect n independent solution vectors x_1(t), ..., x_n(t), and the general solution is their superposition x(t) = c_1 x_1(t) + ... + c_n x_n(t). It is worth bundling these into one object. Stand the n solution vectors up as the columns of a matrix Phi(t) = [x_1(t) | x_2(t) | ... | x_n(t)]. This is the fundamental matrix, and it is the system-level version of the fundamental set of solutions you met in the first guide. With it the general solution collapses to a single matrix-times-vector product, x(t) = Phi(t) c, where c is the column of arbitrary constants.
Because each column solves the system, the whole matrix obeys the system too: dPhi/dt = A Phi, now a matrix differential equation. The test that the columns are genuinely independent is, fittingly, a determinant — det Phi(t), which here plays exactly the role the Wronskian played for single equations. If det Phi is nonzero at one time it is nonzero at every time, and the columns form a basis. This is the same Abel-style all-or-nothing behaviour you saw before, lifted to vectors. A nonsingular fundamental matrix is your certificate that the n trajectories truly span the whole solution space.
The matrix exponential e^{At}
Now for the idea that crowns the subject. The scalar equation x' = a x has the clean closed form x(t) = e^{at} x_0. Could the system dx/dt = A x have a closed form just as clean, x(t) = e^{At} x_0, with the matrix A sitting where the number a sat? The only obstacle is meaning: what could the exponential of a matrix possibly be? We are not raising e to a matrix power literally. Instead we borrow the one definition of the exponential that never mentions roots or repeated multiplication of e — its Taylor series e^z = 1 + z + z^2/2! + z^3/3! + ... — and feed a matrix in. Defining the matrix exponential by e^{At} = I + At + (At)^2/2! + (At)^3/3! + ... makes complete sense: every term is a matrix, you add matrices, and this series provably converges for any square A.
And it does everything you would hope. Differentiate the series term by term in t — legitimate inside its radius of convergence, exactly as for ordinary power series — and you get d/dt[e^{At}] = A + A^2 t + A^3 t^2/2! + ... = A e^{At}. So x(t) = e^{At} x_0 satisfies dx/dt = A x, and at t = 0 it equals e^{0} x_0 = I x_0 = x_0. The matrix exponential is therefore the exact solution of the initial value problem for any starting point x_0 in one stroke. It is also the special fundamental matrix normalized to the identity at t = 0: e^{At} = Phi(t) Phi(0)^{-1}. One object propagates every initial condition forward in time.
The phase plane: reading the flow at a glance
For a 2x2 system the most beautiful payoff is geometric. Drop the time axis and draw the trajectories directly in the (x_1, x_2) plane — the phase plane — each curve the path a solution sweeps out, with little arrows showing the direction of flow. Because A x is constant in time (the system is autonomous), the picture is a fixed vector field, and the origin, where A x = 0, is the lone equilibrium. The astonishing fact is that the two eigenvalues of one 2x2 matrix completely determine the shape of this whole portrait. Their signs and their realness sort every linear flow into a short list of named pictures — the phase-plane classification.
Walk through the cases by feel. Two real eigenvalues of the same sign give a node: if both are negative, every trajectory flows inward and the origin is a stable sink, all paths draining to the center; if both positive, a source, everything flying out. Two real eigenvalues of opposite signs give a saddle — the case in our worked 2x2 above, lambda = 3 and lambda = -1: trajectories rush in along the negative-eigenvalue eigenvector, swerve, and shoot out along the positive one, tracing hyperbola-like paths past an unstable junction. Complex eigenvalues lambda = alpha +/- i beta instead bring rotation, because of the cos and sin hiding inside e^{(alpha + i beta)t}: the real part alpha sets the radius and the imaginary part beta sets the spinning. If alpha is not zero you get a spiral — trajectories wind around the origin while shrinking (alpha < 0, a stable spiral) or growing (alpha > 0, an unstable one), the universal portrait of a damped or driven oscillation winding down or up. If alpha is exactly zero, pure imaginary eigenvalues, you get a center: closed elliptical orbits circling the origin forever, neither growing nor decaying — undamped oscillation, energy perfectly conserved, the phase-plane face of a frictionless pendulum or an ideal LC circuit. Throughout, the sign of the eigenvalue's real part is the verdict on stability: negative real parts pull you home, positive real parts push you away.
- Compute the two eigenvalues from det(A - lambda I) = 0 — equivalently from the trace tau = a + d and determinant Delta = ad - bc of A = [a, b; c, d], since lambda solves lambda^2 - tau lambda + Delta = 0.
- Real and same sign -> node (sink if both negative, source if both positive). Real and opposite signs -> saddle (always unstable).
- Complex with nonzero real part -> spiral (stable if real part < 0, unstable if > 0). Purely imaginary -> center (closed orbits, neutrally stable).
- Read stability off the real parts: both real parts negative means every trajectory returns to the origin; any positive real part means it eventually flees.
Hold the honest boundaries in view. This crisp gallery is exact only for linear systems with constant coefficients; a real nonlinear system has its own, possibly wild, geometry. The bridge back is the linearization you will meet later: near an equilibrium of a nonlinear field you approximate it by its partial-derivative matrix — the Jacobian, an idea straight from multivariable calculus — and the eigenvalues of that local matrix usually predict the nearby portrait. Usually, not always: the borderline cases (a center, a repeated eigenvalue, a zero eigenvalue) are delicate, and there the nonlinear terms you discarded can decide the true behaviour. Within those honest limits, though, you now hold a remarkable power — read two eigenvalues off a matrix, and you can sketch the entire fate of the system without solving a single trajectory.