Differential Equations

A differential equation is an equation that relates a function to its own derivatives. Where an ordinary algebraic equation like x^2 = 9 asks you to find a number, a differential equation asks you to find a whole function. The equation describes how a quantity changes, and solving it recovers the quantity itself. Because rates of change are everywhere in science, differential equations are the natural language of physics, biology, chemistry, engineering, and economics.

A first example

Suppose a population y(t) grows at a rate proportional to its current size. "Rate of change" is the derivative dy/dt, and "proportional to its size" means it equals a constant times y. That sentence, translated, is the differential equation:

dy/dt = k y

A solution is any function that makes the equation true. You will soon verify that y(t) = C e^(k t) works for every constant C. Notice there is not one solution but a whole family, one for each value of C. That family is the general solution. Pinning down C requires one extra fact, an initial condition such as y(0) = 100, which selects the single particular solution that starts at 100.

Order

The order of a differential equation is the highest derivative that appears in it.

• dy/dx = 3x is first order (only a first derivative).
• y'' + 4y = 0 is second order (a second derivative appears).
• y''' - y' = e^x is third order.

This course focuses on first- and second-order equations, which cover an enormous range of applications.

Ordinary versus partial

An ordinary differential equation (ODE) involves a function of a single variable, so all derivatives are ordinary derivatives, like dy/dx. A partial differential equation (PDE) involves a function of several variables and its partial derivatives, like the heat equation. Everything in MATH 308 is an ODE; PDEs are a later course.

Linear versus nonlinear (a preview)

An ODE is linear if the unknown function and its derivatives appear only to the first power and are never multiplied together, with no functions like sin(y) or y^2 wrapped around them. So y'' + 3y' + 2y = cos(x) is linear, but y' = y^2 and y'' + sin(y) = 0 are nonlinear. Linear equations are far more tractable, and most of our exact methods apply to them. We sharpen this definition in the next lesson.

The workflow of the whole course is now visible: read a situation, translate it into a differential equation, solve for the general solution, then use initial conditions to find the particular solution that answers the real question.

Before learning to find solutions, it is worth learning to check them, because checking is always available and never lies. To verify a proposed solution, differentiate it as many times as the equation needs, substitute into both sides, and confirm they are equal for all values of the variable.

Worked example, verifying a solution

Show that y = C e^(3x) solves y' = 3y.

1. Differentiate the candidate: y' = 3 C e^(3x).
2. Compute the right side: 3y = 3 (C e^(3x)) = 3 C e^(3x).
3. The two sides match for every x and every C, so y = C e^(3x) is a solution.

Worked example, a second-order check

Show that y = cos(2x) solves y'' + 4y = 0.

1. First derivative: y' = -2 sin(2x).
2. Second derivative: y'' = -4 cos(2x).
3. Substitute: y'' + 4y = -4 cos(2x) + 4 cos(2x) = 0. True, so it is a solution.

The test for linearity

A differential equation is linear if it can be written in the form

a_n(x) y^(n) + ... + a_1(x) y' + a_0(x) y = g(x)

where the unknown y and its derivatives appear only to the first power, are not multiplied by each other, and sit inside no nonlinear function. The coefficients a_i(x) may depend on x in any way; that does not spoil linearity. A quick checklist for nonlinearity: look for a power like y^2, a product like y y', or a wrapper like sin(y) or e^y. Any one of these makes the equation nonlinear.

Equation	Type	Reason
y' + x^2 y = cos(x)	Linear	y and y' are to the first power; x^2 is just a coefficient
y' = y^2	Nonlinear	y appears squared
y y' = 1	Nonlinear	y is multiplied by y'
y'' + sin(y) = 0	Nonlinear	y sits inside a sine

Initial value problems

An initial value problem (IVP) is a differential equation together with enough initial conditions to fix every arbitrary constant. A first-order equation needs one condition; a second-order equation needs two, typically the value and the first derivative at a starting point.

Worked example. Solve the IVP y' = 3y, y(0) = 5. From the check above the general solution is y = C e^(3x). Apply the condition: y(0) = C e^(0) = C = 5. So the particular solution is y = 5 e^(3x). The initial condition converted an infinite family into the one curve that actually starts at 5.

A first-order equation is separable if you can algebraically move all the y terms (including dy) to one side and all the x terms (including dx) to the other. In symbols, it can be written dy/dx = f(x) g(y), a product of a function of x alone and a function of y alone. Once separated, you integrate each side with respect to its own variable.

The method in three steps

1. Separate: write the equation as dy/g(y) = f(x) dx, gathering y on the left and x on the right.
2. Integrate both sides. A single constant C on one side captures the two constants of integration.
3. Solve for y if possible, then apply any initial condition to find C.

Worked example one

Solve dy/dx = x y.

1. Separate: dy/y = x dx.
2. Integrate: ln|y| = x^2/2 + C.
3. Exponentiate both sides: |y| = e^(x^2/2 + C) = e^C e^(x^2/2). Absorbing +/- e^C into a single constant A gives y = A e^(x^2/2).

So the general solution is y = A e^(x^2/2).

Worked example two, with an initial condition

Solve the IVP dy/dx = -2x y^2, y(0) = 1.

1. Separate: dy/y^2 = -2x dx, that is y^(-2) dy = -2x dx.
2. Integrate: -1/y = -x^2 + C.
3. Solve for y: multiply by -1 to get 1/y = x^2 - C. Writing -C as a new constant D, 1/y = x^2 + D, so y = 1/(x^2 + D).
4. Apply y(0) = 1: 1 = 1/(0 + D), so D = 1.

The particular solution is y = 1/(x^2 + 1). You can check it: y' = -2x/(x^2 + 1)^2 = -2x y^2, exactly the right side.

Watch for lost solutions

Dividing by g(y) is only valid where g(y) is not zero. If g(y_0) = 0 for some constant, then the constant function y = y_0 is also a solution, an equilibrium solution, which the division step can hide. For dy/dx = x y, dividing by y quietly drops the solution y = 0. Always check separately whether any constant that zeroes g(y) is a valid solution.

Many first-order equations are linear but not separable, such as y' + 2y = e^x. For these there is a reliable, universal method: the integrating factor. The whole idea is to multiply the equation by a cleverly chosen function that turns the left side into the derivative of a single product, which you can then undo by integrating.

Step 1, standard form

Write the equation as

y' + P(x) y = Q(x)

with a coefficient of exactly 1 on y'. If some other factor multiplies y', divide the whole equation by it first.

Step 2, the integrating factor

Define the integrating factor

mu(x) = e^(integral of P(x) dx)

The magic of this choice is that after multiplying the equation by mu, the entire left side collapses into d/dx [mu(x) y]. That happens because mu' = P mu by construction, so the product rule reproduces exactly mu y' + P mu y, which is the left side scaled by mu.

Step 3, integrate and solve

The equation becomes d/dx [mu y] = mu Q. Integrate both sides and divide by mu:

y = (1/mu) [ integral of mu(x) Q(x) dx + C ]

Worked example

Solve y' + 2y = e^x.

1. Already in standard form with P(x) = 2 and Q(x) = e^x.
2. Integrating factor: mu = e^(integral of 2 dx) = e^(2x).
3. Multiply through: e^(2x) y' + 2 e^(2x) y = e^(2x) e^x = e^(3x). The left side is d/dx [e^(2x) y].
4. Integrate: e^(2x) y = integral of e^(3x) dx = (1/3) e^(3x) + C.
5. Divide by e^(2x): y = (1/3) e^x + C e^(-2x).

So y = (1/3) e^x + C e^(-2x). The first piece is driven by Q(x); the second, C e^(-2x), is the decaying solution of the associated equation y' + 2y = 0.

Worked example with an initial condition

Solve the IVP x y' + y = 4x, y(1) = 3, for x > 0.

1. Divide by x to reach standard form: y' + (1/x) y = 4. So P = 1/x, Q = 4.
2. Integrating factor: mu = e^(integral of (1/x) dx) = e^(ln x) = x.
3. Multiply: x y' + y = 4x. The left side is d/dx [x y]. (Here the original equation was already in this convenient form.)
4. Integrate: x y = integral of 4x dx = 2x^2 + C, so y = 2x + C/x.
5. Apply y(1) = 3: 2(1) + C = 3, so C = 1.

The particular solution is y = 2x + 1/x. The integrating factor turns any linear first-order equation into a single integration, which is why it is the standard tool.

The simplest and most important model in the whole course says that a quantity changes at a rate proportional to its current amount. Whether it is bacteria dividing, money earning continuous interest, or a radioactive sample decaying, the equation is the same:

dy/dt = k y

This is separable. Separating and integrating, dy/y = k dt gives ln|y| = k t + C, and exponentiating gives the solution

y(t) = y_0 e^(k t)

where y_0 = y(0) is the starting amount. The constant k is the relative growth rate. If k > 0 the quantity grows; if k < 0 it decays.

Doubling time and half-life

For growth, the doubling time is the time for the amount to double. Set y = 2 y_0: 2 y_0 = y_0 e^(k t), so e^(k t) = 2, giving t = ln(2)/k. For decay, the half-life solves e^(k t) = 1/2, giving t = ln(2)/(-k) (positive because k < 0). Both depend only on k, not on the starting amount.

Worked example, population growth

A bacteria culture starts at 500 cells and grows to 800 in 2 hours. Find the model and the population at 5 hours.

1. Model: y(t) = 500 e^(k t).
2. Use the second data point: 800 = 500 e^(2k), so e^(2k) = 1.6, giving 2k = ln(1.6) and k = ln(1.6)/2 approximately 0.235 per hour.
3. At t = 5: y(5) = 500 e^(0.235 * 5) = 500 e^(1.175) approximately 500 (3.238) approximately 1619 cells.

Worked example, radioactive decay

Carbon-14 has a half-life of about 5730 years. Find k and the fraction remaining after 2000 years.

1. Half-life: 1/2 = e^(k * 5730), so k = ln(1/2)/5730 = -ln(2)/5730 approximately -0.000121 per year.
2. After 2000 years: y/y_0 = e^(-0.000121 * 2000) = e^(-0.242) approximately 0.785, so about 78.5 percent remains.

This exact calculation, run in reverse, is the basis of radiocarbon dating. The same equation, one constant, and two data points solve a remarkable range of real problems.

Two classic models push our first-order methods a little further. Both lead to linear equations that also happen to be separable.

Newton's law of cooling

Newton's law of cooling says an object's temperature changes at a rate proportional to the difference between the object and its surroundings. If T(t) is the object's temperature and T_s the constant surrounding temperature, then

dT/dt = -k (T - T_s), with k > 0.

The minus sign makes a hot object cool and a cold object warm, always toward T_s. Separating variables with u = T - T_s (so du = dT) gives du/u = -k dt, hence u = u_0 e^(-k t), that is

T(t) = T_s + (T_0 - T_s) e^(-k t).

As t grows, e^(-k t) to 0, so T to T_s: the object approaches room temperature, exactly as experience predicts.

Worked example. A cup of coffee at 90 degrees sits in a 20 degree room and cools to 70 degrees in 10 minutes. Find its temperature after 20 minutes.

1. Model: T = 20 + (90 - 20) e^(-k t) = 20 + 70 e^(-k t).
2. Use T(10) = 70: 70 = 20 + 70 e^(-10k), so 50 = 70 e^(-10k), giving e^(-10k) = 5/7 and k = -ln(5/7)/10 = ln(7/5)/10 approximately 0.0336.
3. At t = 20: e^(-20k) = (e^(-10k))^2 = (5/7)^2 = 25/49, so T(20) = 20 + 70 (25/49) approximately 20 + 35.7 = 55.7 degrees.

Mixing problems

Imagine a tank holding a well-stirred solution, with liquid flowing in and out. Let A(t) be the amount of solute (say salt) in the tank. The governing principle is a balance:

dA/dt = (rate in) - (rate out).

Each rate is a concentration times a flow rate. The rate in is the inflow concentration times the inflow speed. The rate out is the tank's current concentration, A/V where V is the volume, times the outflow speed.

Worked example. A 100 liter tank holds pure water. Brine with 2 grams of salt per liter flows in at 5 liters per minute, and the well-mixed solution flows out at 5 liters per minute (so volume stays 100 L). Find A(t).

1. Rate in: 2 grams/L * 5 L/min = 10 grams/min.
2. Rate out: concentration A/100 times 5, that is 5A/100 = A/20 grams/min.
3. Equation: dA/dt = 10 - A/20. In standard linear form, A' + (1/20) A = 10.
4. Integrating factor mu = e^(t/20): d/dt [e^(t/20) A] = 10 e^(t/20), so e^(t/20) A = 200 e^(t/20) + C, giving A = 200 + C e^(-t/20).
5. Initial condition A(0) = 0 (pure water): 0 = 200 + C, so C = -200.

Thus A(t) = 200 (1 - e^(-t/20)). As t to infinity, A to 200 grams, which matches intuition: the tank eventually reaches the inflow concentration of 2 grams per liter times 100 liters.

Second-order linear equations model any system with inertia and a restoring force, from a vibrating spring to an electrical circuit. The general constant-coefficient form is

a y'' + b y' + c y = g(x)

with constants a, b, c and a not zero. The equation is homogeneous when g(x) = 0, meaning there is no external driving term, and nonhomogeneous when g(x) is present. This module handles the homogeneous case, the essential foundation.

Why two solutions, combined

A second-order linear homogeneous equation has a beautiful structure. If y_1 and y_2 are both solutions, then so is every combination C_1 y_1 + C_2 y_2. This is the superposition principle, and it holds precisely because the equation is linear and homogeneous. To capture all solutions we need two building-block solutions that are genuinely different, not multiples of each other.

Two functions are linearly independent if neither is a constant multiple of the other. For instance e^x and e^(2x) are independent, but e^x and 3 e^x are not. The general solution of a second-order homogeneous equation is

y = C_1 y_1 + C_2 y_2

where y_1 and y_2 are any two linearly independent solutions. Two arbitrary constants is exactly right for a second-order equation, and they are the two facts an initial value problem supplies, usually y(0) and y'(0).

The characteristic equation

How do we find y_1 and y_2? Guess that a solution has the form y = e^(r x) for some number r. Then y' = r e^(r x) and y'' = r^2 e^(r x). Substitute into a y'' + b y' + c y = 0:

a r^2 e^(r x) + b r e^(r x) + c e^(r x) = 0.

Since e^(r x) is never zero, divide it out to leave a plain quadratic in r, the characteristic equation:

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

The exponents in the solution are just the roots of this quadratic. Everything now depends on what kind of roots the quadratic has, which the discriminant b^2 - 4ac decides. The next lesson works through all three cases. This translation, from a differential equation to an ordinary quadratic, is the central trick of the whole module.

The general solution of a y'' + b y' + c y = 0 is decided entirely by the roots of its characteristic equation a r^2 + b r + c = 0. The discriminant b^2 - 4ac sorts every problem into one of three cases.

Case 1, two distinct real roots

If b^2 - 4ac > 0, the quadratic has two different real roots r_1 and r_2. Each gives an independent exponential solution, so

y = C_1 e^(r_1 x) + C_2 e^(r_2 x).

Worked example. Solve y'' - 5y' + 6y = 0. The characteristic equation r^2 - 5r + 6 = 0 factors as (r - 2)(r - 3) = 0, so r = 2, 3. The general solution is y = C_1 e^(2x) + C_2 e^(3x).

Case 2, one repeated real root

If b^2 - 4ac = 0, the quadratic has a single repeated root r = -b/(2a). That gives only one exponential, but a second-order equation needs two independent solutions. The missing partner is the same exponential multiplied by x. So

y = C_1 e^(r x) + C_2 x e^(r x).

Worked example. Solve y'' - 6y' + 9y = 0. Here r^2 - 6r + 9 = (r - 3)^2 = 0, a repeated root r = 3. The general solution is y = C_1 e^(3x) + C_2 x e^(3x). You can verify x e^(3x) is a genuine second solution by substitution.

Case 3, complex conjugate roots

If b^2 - 4ac < 0, the roots are a complex conjugate pair r = alpha +/- beta i, where alpha = -b/(2a) and beta = sqrt(4ac - b^2)/(2a). Using Euler's formula to convert the complex exponentials into real trig functions, the real general solution is

y = e^(alpha x) [C_1 cos(beta x) + C_2 sin(beta x)].

The real part alpha controls growth or decay of the envelope, and beta is the frequency of oscillation. This is the case that produces vibrations and waves.

Worked example. Solve y'' + 4y = 0. The characteristic equation r^2 + 4 = 0 gives r = +/- 2i, so alpha = 0 and beta = 2. The general solution is y = C_1 cos(2x) + C_2 sin(2x), pure oscillation with no decay.

Worked example with decay. Solve y'' + 2y' + 5y = 0. The characteristic equation r^2 + 2r + 5 = 0 has discriminant 4 - 20 = -16 < 0, and roots r = (-2 +/- sqrt(-16))/2 = -1 +/- 2i. So alpha = -1, beta = 2, and y = e^(-x)[C_1 cos(2x) + C_2 sin(2x)], a decaying oscillation.

Discriminant	Roots	General solution
b^2 - 4ac > 0	real, distinct r1, r2	C1 e^(r1 x) + C2 e^(r2 x)
b^2 - 4ac = 0	repeated r	C1 e^(r x) + C2 x e^(r x)
b^2 - 4ac < 0	alpha +/- beta i	e^(alpha x)[C1 cos(beta x) + C2 sin(beta x)]

Now restore the driving term: a y'' + b y' + c y = g(x) with g(x) not zero. The key structural fact is that the general solution splits into two parts:

y = y_h + y_p

where y_h is the homogeneous solution (the general solution of the equation with g = 0, from Module 4) and y_p is any single particular solution of the full equation. The homogeneous part carries the two arbitrary constants; the particular part answers the forcing.

The method of undetermined coefficients

When g(x) is a polynomial, an exponential, a sine or cosine, or a product of these, we can guess the form of y_p with unknown coefficients, then solve for them. The trial form mirrors the driving term:

Driving term g(x)	Trial y_p
polynomial, e.g. 3x + 2	A x + B (same degree, all terms)
exponential, e.g. e^(2x)	A e^(2x)
sine or cosine, e.g. cos(3x)	A cos(3x) + B sin(3x)
product, e.g. x e^x	(A x + B) e^x

Worked example, exponential forcing

Solve y'' - 5y' + 6y = e^x.

1. Homogeneous part: from Module 4, y_h = C_1 e^(2x) + C_2 e^(3x).
2. Trial: g = e^x, so guess y_p = A e^x. Then y_p' = A e^x and y_p'' = A e^x.
3. Substitute: A e^x - 5 A e^x + 6 A e^x = e^x, so (1 - 5 + 6) A e^x = 2 A e^x = e^x.
4. Solve: 2A = 1, so A = 1/2 and y_p = (1/2) e^x.

The full solution is y = C_1 e^(2x) + C_2 e^(3x) + (1/2) e^x.

Worked example, trig forcing

Solve y'' + y = 4 cos(x)... but first a warning. The homogeneous solution is y_h = C_1 cos(x) + C_2 sin(x), which already contains cos(x). When the trial form duplicates a homogeneous solution, multiply the trial by x. So use y_p = x (A cos x + B sin x). Differentiating and substituting (a longer computation) yields A = 0, B = 2, so y_p = 2x sin(x) and y = C_1 cos x + C_2 sin x + 2x sin x. This growing term is resonance, which we meet again with springs.

A cleaner trig example

Solve y'' - 5y' + 6y = 3 cos(x). No overlap with y_h, so try y_p = A cos x + B sin x. Then y_p' = -A sin x + B cos x and y_p'' = -A cos x - B sin x. Substituting and grouping the cos x and sin x terms gives the system 5A - 5B = 3 and 5A + 5B = 0. Adding: 10A = 3, so A = 3/10; subtracting: B = -3/10. Thus y_p = (3/10) cos x - (3/10) sin x, and the full solution is y = C_1 e^(2x) + C_2 e^(3x) + (3/10) cos x - (3/10) sin x.

The reason second-order linear equations matter so much is that one equation describes both a bouncing spring and an electrical circuit. Learning the mechanics teaches the circuit for free.

The spring-mass system

A mass m on a spring with stiffness k, with a damping (friction) force proportional to velocity with constant c, and an external force F(t), obeys Newton's second law:

m x'' + c x' + k x = F(t)

Here x(t) is displacement from equilibrium. Each term is a force: m x'' is inertia, c x' is damping, and k x is the spring's restoring pull (Hooke's law). This is exactly the constant-coefficient form of Module 4, so the characteristic-equation method applies directly.

Free undamped motion

With no damping and no forcing (c = 0, F = 0), m x'' + k x = 0 gives roots r = +/- i sqrt(k/m), so

x(t) = C_1 cos(omega_0 t) + C_2 sin(omega_0 t), with natural frequency omega_0 = sqrt(k/m).

The mass oscillates forever at its natural frequency. Worked example: m = 1, k = 16 gives omega_0 = sqrt(16) = 4, so x = C_1 cos(4t) + C_2 sin(4t).

Damping, three regimes

Turning on damping (c > 0) determines behavior through the discriminant c^2 - 4mk:

• Underdamped (c^2 < 4mk): complex roots, so the mass oscillates with an amplitude that decays like e^(-(c/2m) t). A struck tuning fork.
• Critically damped (c^2 = 4mk): a repeated root, the fastest return to rest with no oscillation. A well-designed door closer.
• Overdamped (c^2 > 4mk): two negative real roots, a slow sluggish return to rest. Moving through thick oil.

Worked example. x'' + 2x' + 5x = 0 models m = 1, c = 2, k = 5. The discriminant is 4 - 20 = -16 < 0, underdamped, with roots -1 +/- 2i, so x = e^(-t)[C_1 cos(2t) + C_2 sin(2t)]: oscillation shrinking under an e^(-t) envelope.

Resonance

If an undamped system is driven at its natural frequency, the particular solution grows without bound (the x-times factor from the last lesson). This is resonance, the reason soldiers break step on bridges and a pushed swing climbs higher with each well-timed push.

The RLC circuit, same equation

A series circuit with inductance L, resistance R, and capacitance C, driven by a voltage E(t), satisfies, for the charge q(t) on the capacitor,

L q'' + R q' + (1/C) q = E(t).

Compare term by term: inductance L plays the role of mass (inertia of current), resistance R plays the role of damping, and 1/C plays the role of spring stiffness. An underdamped circuit rings; an overdamped one settles quietly. One mathematical solution, two physical worlds.

The Laplace transform is a machine that converts a function of time f(t) into a new function of a variable s. Its power is that calculus operations (derivatives) become algebra operations (multiplication) in the new world, so a differential equation becomes an algebraic equation you can solve by hand. The transform is defined by the integral

F(s) = L{f(t)} = integral from 0 to infinity of e^(-s t) f(t) dt

for values of s large enough that the integral converges. We write L{f(t)} = F(s), using a capital letter for the transform.

Worked example, from the definition

Find L{1}. With f(t) = 1:

1. L{1} = integral from 0 to infinity of e^(-s t) dt.
2. The antiderivative is -(1/s) e^(-s t). Evaluate from 0 to infinity: at the top e^(-s t) to 0 (for s > 0), at the bottom it is -(1/s)(1) = -1/s.
3. The result is 0 - (-1/s) = 1/s.

So L{1} = 1/s. A similar computation gives L{e^(a t)} = 1/(s - a), valid for s > a.

A table of essential transforms

These are the building blocks; memorize the common ones and keep the rest handy.

f(t)	F(s) = L{f(t)}
1	1/s
t	1/s^2
t^n	n! / s^(n+1)
e^(a t)	1/(s - a)
sin(b t)	b/(s^2 + b^2)
cos(b t)	s/(s^2 + b^2)

Linearity

The Laplace transform is linear: the transform of a sum is the sum of the transforms, and constants pull out.

L{a f(t) + b g(t)} = a L{f(t)} + b L{g(t)}.

Worked example. Find L{3 + 2t - e^(4t)}. Apply the table term by term: L{3} = 3/s, L{2t} = 2/s^2, and L{e^(4t)} = 1/(s - 4). So the transform is

3/s + 2/s^2 - 1/(s - 4).

Linearity plus the short table already lets you transform any polynomial, exponential, or sinusoid. In the next lesson we transform derivatives, which is where the method earns its keep.

The Laplace method shines on initial value problems because the transform of a derivative automatically absorbs the initial conditions. The two key rules are

L{y'} = s Y(s) - y(0)

L{y''} = s^2 Y(s) - s y(0) - y'(0)

where Y(s) = L{y(t)}. Each derivative becomes multiplication by s, minus the initial data. A differential equation in t turns into an algebraic equation in Y(s).

The three-step recipe

1. Transform both sides of the equation, using the derivative rules to insert the initial conditions.
2. Solve algebraically for Y(s).
3. Invert Y(s) back to y(t) using the table, with partial fractions if needed.

Worked example, first order

Solve y' + 3y = 0, y(0) = 2.

1. Transform: [s Y - y(0)] + 3 Y = 0, that is s Y - 2 + 3 Y = 0.
2. Solve: (s + 3) Y = 2, so Y = 2/(s + 3).
3. Invert: 1/(s + 3) is the transform of e^(-3t), so y = 2 e^(-3t).

Worked example, second order

Solve y'' + 4y = 0, y(0) = 3, y'(0) = 0.

1. Transform: [s^2 Y - s y(0) - y'(0)] + 4 Y = 0, that is s^2 Y - 3s - 0 + 4 Y = 0.
2. Solve: (s^2 + 4) Y = 3s, so Y = 3s/(s^2 + 4).
3. Invert: s/(s^2 + 4) is the transform of cos(2t), so y = 3 cos(2t).

Notice the initial conditions were built in from the start, so there is no separate step of solving for arbitrary constants at the end. That is the practical advantage over the characteristic-equation method.

Partial fractions for the inverse

Often Y(s) is a rational function that is not directly in the table, and partial fractions break it into pieces that are.

Worked example. Solve y' - y = e^(2t), y(0) = 0.

1. Transform: s Y - 0 - Y = 1/(s - 2), so (s - 1) Y = 1/(s - 2) and Y = 1/[(s - 1)(s - 2)].
2. Partial fractions: 1/[(s - 1)(s - 2)] = A/(s - 1) + B/(s - 2). Clearing denominators, 1 = A(s - 2) + B(s - 1). Set s = 2: 1 = B. Set s = 1: 1 = -A, so A = -1.
3. So Y = -1/(s - 1) + 1/(s - 2), which inverts term by term to y = -e^t + e^(2t).

You can check: y(0) = -1 + 1 = 0 matches, and y' - y = (-e^t + 2e^(2t)) - (-e^t + e^(2t)) = e^(2t), the correct right side.

Real systems often have several quantities changing together, each rate depending on all of them. Two interacting populations, or two coupled tanks, lead to a system of first-order equations. A linear system in two unknowns x(t) and y(t) looks like

x' = a x + b y, y' = c x + d y,

which we write compactly in matrix form as X' = A X, where X = [x, y] and A is the 2x2 coefficient matrix with rows [a, b] and [c, d].

The eigenvalue idea

Just as the single equation y' = k y had solution e^(k t), we look for solutions of the form X = v e^(lambda t), where v is a constant vector and lambda a constant. Substituting gives lambda v e^(lambda t) = A v e^(lambda t), and cancelling the exponential leaves the eigenvalue equation

A v = lambda v.

So lambda must be an eigenvalue of A and v a matching eigenvector. Eigenvalues come from the characteristic equation det(A - lambda I) = 0, which for a 2x2 matrix is

lambda^2 - (a + d) lambda + (a d - b c) = 0.

Worked example

Solve the system x' = x + y, y' = 4x + y. The matrix A has rows [1, 1] and [4, 1].

1. Eigenvalues: det(A - lambda I) = (1 - lambda)^2 - (1)(4) = 0. Expand: (1 - lambda)^2 = 4, so 1 - lambda = +/- 2, giving lambda = -1 and lambda = 3.
2. Eigenvector for lambda = 3: solve (A - 3I) v = 0. The first row gives (1 - 3) v_1 + v_2 = -2 v_1 + v_2 = 0, so v_2 = 2 v_1; take v = [1, 2].
3. Eigenvector for lambda = -1: the first row gives (1 + 1) v_1 + v_2 = 2 v_1 + v_2 = 0, so v_2 = -2 v_1; take v = [1, -2].

The general solution combines the two eigen-solutions:

X = C_1 [1, 2] e^(3t) + C_2 [1, -2] e^(-t).

Written out, x = C_1 e^(3t) + C_2 e^(-t) and y = 2 C_1 e^(3t) - 2 C_2 e^(-t). The sign of each eigenvalue tells the story: the lambda = 3 direction grows, the lambda = -1 direction decays. Real distinct eigenvalues behave much like the two-real-roots case of a single second-order equation, which is no accident: a second-order equation can always be rewritten as a first-order system.

Most differential equations that arise in practice cannot be solved with a neat formula. When exact methods fail, we approximate the solution numerically, computing the value of y at a sequence of points. The simplest such scheme is Euler's method, and it captures the core idea behind all the others.

The idea, follow the slope

Suppose we have an initial value problem y' = f(x, y) with y(x_0) = y_0. The equation hands us the slope at any point: it is just f(x, y). Euler's method starts at the known point and repeatedly takes a short straight step in the direction of the current slope. Choose a small step size h. Then update:

x_(n+1) = x_n + h

y_(n+1) = y_n + h f(x_n, y_n)

In words: the new y is the old y plus the step size times the slope there. Geometrically, you walk along the tangent line for a short distance, then recompute the slope and repeat, tracing a polygonal path that approximates the true curve.

Worked example

Approximate y(0.3) for y' = x + y, y(0) = 1, using step size h = 0.1. Here f(x, y) = x + y.

n	x_n	y_n	slope f = x_n + y_n	y_(n+1) = y_n + 0.1 f
0	0.0	1.000	0 + 1.000 = 1.000	1.000 + 0.1(1.000) = 1.100
1	0.1	1.100	0.1 + 1.100 = 1.200	1.100 + 0.1(1.200) = 1.220
2	0.2	1.220	0.2 + 1.220 = 1.420	1.220 + 0.1(1.420) = 1.362
3	0.3	1.362	-	-

So Euler's method estimates y(0.3) approximately 1.362. The exact solution of this equation is y = 2 e^x - x - 1, which gives y(0.3) = 2 e^(0.3) - 1.3 approximately 2(1.34986) - 1.3 approximately 1.400. Our estimate 1.362 is close but low, an error of about 0.038.

Step size and error

Euler's method is only approximate because it assumes the slope stays constant across each whole step, when in reality the slope changes continuously. Two consequences follow:

• Smaller h is more accurate. Halving the step size roughly halves the error of Euler's method, at the cost of twice as many steps. Euler's method is called first-order accurate for this reason.
• Error accumulates. Each step's small error carries into the next, so estimates drift over long intervals.

More refined schemes, such as the improved Euler (Heun) method and the widely used Runge-Kutta methods, reduce error dramatically by sampling the slope at several points within each step before committing. But every one of them rests on the same foundation you just used: read the slope from the equation, take a controlled step, and repeat. That single idea is how computers solve the differential equations that no formula can.