Higher-Order Differential Equations

Second-order and higher differential equations look intimidating, but the MSTE paper almost always hands you the friendly case: linear, with constant coefficients. Every one of those reduces to a two-step recipe. First solve the characteristic equation for its roots and read off the homogeneous solution; then, if a forcing term is present, add a particular solution built by undetermined coefficients. This lesson walks the whole pipeline, from writing the characteristic equation to fixing the arbitrary constants with initial or boundary conditions, and carries a damped-oscillation problem and a polynomial-forcing problem all the way to a numeric answer.

Linear constant-coefficient ODEs and the characteristic equation

A linear ODE with constant coefficients has the shape a y'' + b y' + c y = g(x), where a, b, and c are numbers, not functions of x. When g(x) = 0 the equation is homogeneous; when g(x) is nonzero it is nonhomogeneous and g(x) is the forcing term.

The trick that unlocks the homogeneous case is to guess y = e^(rx). Because every derivative of e^(rx) just multiplies it by another factor of r, substituting turns the ODE into a plain polynomial in r. For a y'' + b y' + c y = 0 you replace y'' with r^2, y' with r, and y with 1, giving the characteristic equation:

a r^2 + b r + c = 0.

Solve that quadratic (or higher polynomial for higher order) and its roots hand you the building blocks of the solution. A second-order equation always needs two independent building blocks, and in general an n-th order linear equation needs n of them.

Roots of the characteristic equation

The nature of the roots, set by the discriminant b^2 - 4ac, decides the form of the homogeneous solution. There are exactly three cases, and memorizing this table is most of the battle.

Discriminant b^2 - 4ac	Roots	Homogeneous solution y_h
> 0	two real distinct roots r1, r2	y = C1 e^(r1 x) + C2 e^(r2 x)
= 0	one repeated real root r	y = (C1 + C2 x) e^(r x)
< 0	complex pair alpha +/- beta i	y = e^(alpha x) (C1 cos(beta x) + C2 sin(beta x))

The same three cases describe a spring-mass-damper: two real roots mean an overdamped return with no oscillation, a repeated root means the critically damped borderline (the fastest non-oscillating return), and a complex pair means an underdamped, oscillating decay.

Homogeneous solutions from the roots

Once you have the roots, writing y_h is mechanical. Consider y'' - 5y' + 6y = 0. The characteristic equation is r^2 - 5r + 6 = 0, which factors as (r - 2)(r - 3) = 0, so r = 2 and r = 3. These are two real distinct roots, giving y_h = C1 e^(2x) + C2 e^(3x).

Worked example (repeated root): Solve y'' - 6y' + 9y = 0 with y(0) = 4 and y'(0) = 5. The characteristic equation r^2 - 6r + 9 = 0 is a perfect square (r - 3)^2 = 0, so r = 3 is a repeated root. The solution form is y = (C1 + C2 x) e^(3x). Apply y(0): C1 = 4. Differentiate, y' = C2 e^(3x) + 3(C1 + C2 x) e^(3x), so y'(0) = C2 + 3 C1 = 5, which gives C2 = 5 - 12 = -7. The particular solution is y = (4 - 7x) e^(3x), and evaluating at x = 1 gives y(1) = (4 - 7) e^3 = -3 e^3, about -60.3 units.

A repeated root supplies only one function, e^(3x). The second independent piece must be x e^(3x); without that extra factor of x you would be one solution short of the two a second-order equation demands.

Complex roots and damped oscillation

When the discriminant is negative the roots come as a conjugate pair alpha +/- beta i. Euler's formula converts the complex exponentials into real sines and cosines, and the solution becomes y = e^(alpha x)(C1 cos(beta x) + C2 sin(beta x)). Read the two parts physically: the real part alpha is the exponential envelope that shrinks (alpha < 0) or grows (alpha > 0) the amplitude, and the imaginary part beta is the angular frequency of the oscillation inside that envelope.

Worked example (damped oscillation carried to a number): A 1 kg mass on a spring of stiffness k = 20 N/m with a damper of coefficient c = 4 N-s/m obeys m y'' + c y' + k y = 0, that is y'' + 4 y' + 20 y = 0, where y is displacement in metres and t is in seconds. The characteristic equation r^2 + 4r + 20 = 0 has discriminant 4^2 - 4(1)(20) = 16 - 80 = -64, which is negative, so the roots are r = (-4 +/- sqrt(-64)) / 2 = (-4 +/- 8i) / 2 = -2 +/- 4i. Hence alpha = -2 and beta = 4, and c^2 - 4mk = -64 < 0 confirms the system is underdamped. The solution is y(t) = e^(-2t)(C1 cos(4t) + C2 sin(4t)). Release it from y(0) = 0.1 m with y'(0) = 0. From y(0): C1 = 0.1 m. Differentiating, y'(t) = e^(-2t)[(-2C1 + 4C2) cos(4t) + (-2C2 - 4C1) sin(4t)], so y'(0) = -2C1 + 4C2 = 0, giving C2 = C1 / 2 = 0.05 m. Therefore y(t) = e^(-2t)(0.1 cos(4t) + 0.05 sin(4t)) m. At t = 0.5 s, 4t = 2 rad, so cos(2) = -0.4161 and sin(2) = 0.9093, and y(0.5) = e^(-1)(0.1(-0.4161) + 0.05(0.9093)) = 0.3679(0.003850) = 0.00142 m, that is about 1.4 mm.

Method of undetermined coefficients

For a nonhomogeneous equation you still need a particular solution y_p, one function that satisfies the full equation forcing and all. When the forcing g(x) is a polynomial, an exponential, a sine or cosine, or a product of those, undetermined coefficients works: guess a trial y_p that mirrors the shape of g(x), plug it in, and match coefficients to solve for the unknown constants.

Forcing term g(x)	Trial particular y_p
constant k	A
polynomial of degree n	A_n x^n + ... + A_1 x + A_0
e^(rx)	A e^(rx)
sin(bx) or cos(bx)	A cos(bx) + B sin(bx)
e^(rx) sin(bx) or e^(rx) cos(bx)	e^(rx)(A cos(bx) + B sin(bx))

One catch, the modification rule: if any term of your trial y_p already appears in the homogeneous solution y_h, multiply the whole trial by x (or by x^2 for a doubled root) until nothing overlaps. For example if g(x) = e^(2x) and r = 2 is a simple root of the characteristic equation, then A e^(2x) is already part of y_h, so the correct trial is A x e^(2x).

Particular and general solutions with a polynomial forcing

Worked example (polynomial forcing carried to a number): Solve y'' + 3 y' + 2 y = 2 t^2, where y is in metres and t is in seconds, with y(0) = 0 and y'(0) = 0. Start with the homogeneous side: r^2 + 3r + 2 = 0 factors as (r + 1)(r + 2) = 0, so r = -1 and r = -2, and y_h = C1 e^(-t) + C2 e^(-2t). The forcing 2 t^2 is a degree-2 polynomial, so the trial is y_p = A t^2 + B t + C, giving y_p' = 2A t + B and y_p'' = 2A. Substitute: 2A + 3(2A t + B) + 2(A t^2 + B t + C) = 2 t^2. Grouping by power of t, the t^2 terms give 2A = 2 so A = 1; the t^1 terms give 6A + 2B = 0 so B = -3; the constant terms give 2A + 3B + 2C = 0, that is 2 - 9 + 2C = 0, so C = 3.5. Thus y_p = t^2 - 3t + 3.5, and the general solution is y = C1 e^(-t) + C2 e^(-2t) + t^2 - 3t + 3.5.

Now fit the initial conditions to the full general solution, never to y_h alone. From y(0) = 0: C1 + C2 + 3.5 = 0, so C1 + C2 = -3.5. Differentiate, y' = -C1 e^(-t) - 2C2 e^(-2t) + 2t - 3, and y'(0) = 0 gives -C1 - 2C2 - 3 = 0, so C1 + 2C2 = -3. Subtracting the two, C2 = 0.5 and C1 = -4. The complete solution is y(t) = -4 e^(-t) + 0.5 e^(-2t) + t^2 - 3t + 3.5 m. Check it at t = 1 s: y(1) = -4 e^(-1) + 0.5 e^(-2) + 1 - 3 + 3.5 = -4(0.3679) + 0.5(0.1353) + 1.5 = -1.4716 + 0.0677 + 1.5 = 0.096 m, that is about 96 mm. The two exponentials are the transient that fades as t grows, and t^2 - 3t + 3.5 is the lasting particular part the solution settles onto.

Initial versus boundary conditions

The general solution of an n-th order equation carries n arbitrary constants, and you need n conditions to nail them down. How those conditions are placed decides the label. An initial-value problem (IVP) gives all conditions at a single point, typically y and y' at t = 0, which is the natural setup for a system starting from a known state and running forward in time. A boundary-value problem (BVP) instead gives conditions at two different points, for instance y(0) = 0 and y(L) = 0 for a beam pinned at both ends. Solve the general solution first in both cases, then substitute the conditions to build and solve a small system for the constants. Apply the conditions last, after y_h and y_p are already combined into the full general solution.

Exam-day strategy

• Confirm the equation is linear with constant coefficients before reaching for the characteristic equation; a variable coefficient or a y^2 term means this method does not apply.
• Write the characteristic equation by swapping y'' for r^2, y' for r, and y for 1, then factor or use the quadratic formula; check the discriminant sign to pick the real, repeated, or complex case.
• For a repeated root do not write C1 e^(rx) + C2 e^(rx); that is one solution twice. The second piece must carry a factor of x, giving (C1 + C2 x) e^(rx).
• On complex roots alpha +/- beta i, read alpha as the decay or growth in the exponential envelope and beta as the oscillation frequency, then write e^(alpha x)(C1 cos(beta x) + C2 sin(beta x)).
• Match the trial y_p to the shape of the forcing g(x); if any trial term already lives in y_h, multiply the trial by x (or x^2 for a doubled root) before solving for the coefficients.
• Assemble the full general solution y = y_h + y_p first, then apply the initial or boundary conditions to solve for the constants; fitting the conditions to y_h alone is the single most common wrong turn.