Chapter 8
Ordinary Differential Equations
8.1Introduction and General Form
An ordinary differential equation (ODE) is an equation that relates a function \(y(x)\) to its derivatives with respect to a single independent variable \(x\):
Order of an ODE: The order of a differential equation is the highest derivative that appears in the equation.
\begin{equation} \frac{d^3y}{dx^3} + 2\frac{dy}{dx} = 0 \quad \Rightarrow \quad \text{Order} = 3 \end{equation}Degree of an ODE: The degree of a differential equation is the exponent of the highest-order derivative, provided the equation is polynomial in derivatives.
\begin{equation} \left(\frac{d^2y}{dx^2}\right)^3 + \frac{dy}{dx} = 0 \quad \Rightarrow \quad \text{Degree} = 3 \end{equation}Solution of an ODE: A function \(y(x)\) that satisfies the differential equation for all \(x\) in a given interval.
\begin{equation} \frac{dy}{dx} = 2x \quad \Rightarrow \quad y = x^2 + C \end{equation}
8.2Classification of ODEs
Ordinary differential equations (ODEs) can be classified according to several characteristics. The most fundamental distinction is between linear and nonlinear equations.
Linear vs. Nonlinear ODEs
An ODE is called linear if the dependent variable \(y\) and all its derivatives appear only to the first power and are not multiplied or composed with each other. The general \(n^{\text{th}}\)-order linear ODE can be written as
where \(a_i(x)\) and \(f(x)\) are known functions of the independent variable \(x\).
If any term involves products or nonlinear functions of \(y\) or its derivatives (such as \(y^2\), \(y\,y'\), or \(\sin y\)), the equation is said to be nonlinear.
Linear ODEs are generally easier to solve because they obey the principle of superposition, if \(y_1\) and \(y_2\) are solutions, then any linear combination \(C_1 y_1 + C_2 y_2\) is also a solution. Nonlinear ODEs, on the other hand, can exhibit far richer behavior, such as limited analytic solvability, sensitivity to initial conditions, or even chaotic motion. However, certain first-order nonlinear equations, such as separable or Bernoulli equations, can still be solved analytically.
Classification of Linear ODEs
Once an equation is identified as linear, we can further classify it based on its right-hand side and coefficients.
Homogeneous Linear ODE: The right-hand side equals zero,
\begin{equation} a_n(x)\,y^{(n)} + a_{n-1}(x)\,y^{(n-1)} + \cdots + a_0(x)\,y = 0. \end{equation}In this case, every term involves \(y\) or its derivatives, and the trivial solution \(y=0\) always satisfies the equation. The general solution typically involves \(n\) arbitrary constants determined by initial or boundary conditions.
Nonhomogeneous Linear ODE: The right-hand side is a nonzero function \(f(x)\),
\begin{equation} a_n(x)\,y^{(n)} + a_{n-1}(x)\,y^{(n-1)} + \cdots + a_0(x)\,y = f(x). \end{equation}The general solution is the sum of two parts:
\begin{equation} \boxed{y(x) = y_h(x) + y_p(x)}, \end{equation}where \(y_h\) is the general solution of the corresponding homogeneous equation and \(y_p\) is a particular solution to the full nonhomogeneous equation.
Linear ODEs with Constant Coefficients: When the coefficients \(a_i\) are constants rather than functions of \(x\),
\begin{equation} a_n\,y^{(n)} + a_{n-1}\,y^{(n-1)} + \cdots + a_1\,y' + a_0\,y = f(x), \end{equation}the equation becomes especially important in physics and engineering. Such forms describe systems like simple harmonic oscillators, RC and RL circuits, and exponential decay processes.
8.3First-Order ODEs
A first-order ODE involves only the first derivative of the unknown function:
This general form can represent both linear and nonlinear cases, depending on how \(f(x, y)\) depends on \(y\). Such equations appear widely in physics—for example, in modeling exponential decay, population growth, or motion with velocity-dependent forces.
8.3.1Separable Equations
A first-order ODE is called separable if it can be expressed as
where the function on the right-hand side is the product of a function of \(x\) and a function of \(y\). Although separable equations are often nonlinear, they are still analytically solvable because the variables can be separated.
Rearranging terms gives
and integrating both sides yields
Solution.
Step 1. Recognize the type of equation. This is a first-order separable (and also linear) differential equation.
Step 2. Separate the variables.
Step 3. Integrate both sides.
Step 4. Simplify the result. Exponentiate both sides:
Let \(C_1 = e^{C}\) (a new arbitrary constant). Then
Step 5. Verification. Differentiate \(y = C_1 e^{x^{2}/2}\):
confirming that the solution satisfies the ODE.
Solution.
Step 1. Separate the variables. Divide both sides by \(x(y+1)\):
Step 2. Integrate both sides.
Step 3. Simplify the constant. Let \(C = \ln a\) (for \(a>0\)) to combine logarithms conveniently:
Step 4. Remove the logarithms.
Because \(a\) can absorb the sign, we write simply
Step 5. Interpretation. The general solution
represents a family of straight lines in the \((x,y)\)-plane, one line for each value of the constant \(a\). Selecting a particular value of \(a\) determines a single specific curve from this family.
Find the particular solutions satisfying (a) \(x=0, y=0\) and (b) \(x=0, y=1\), and discuss any solutions not obtained by separation.
Solution.
Step 1. Separation of variables and integration. For \(|y| < 1\), the equation can be written in separable form:
Integrating both sides gives
so the general solution (valid where \(|y|<1\)) is
Step 2. Monotonicity and physical interpretation. Since \(y' = \sqrt{1 - y^2} \ge 0\), the slope of the curve is always nonnegative. Hence, only the increasing portions of the sine curve correspond to valid solution segments. That is, we must restrict to intervals where
The figure below shows the slope field and the increasing sine segments that satisfy this condition.

Step 3. Singular (boundary) solutions. In separating variables, we divided by \(\sqrt{1 - y^2}\), which is not defined when \(|y| = 1\). However, substituting directly into the differential equation shows that
are valid constant solutions since \(y' = 0 = \sqrt{1 - y^2}\). These are singular solutions that cannot be obtained from the separated form.
Step 4. Particular solutions.
(a) \(y(0)=0\): From \(\arcsin y = x + \alpha\), substituting \(x=0,\, y=0\) gives \(\alpha=0\) or \(\pi\). Because \(y'\ge 0\), we take the increasing branch:
To extend this continuously for all \(x\), attach the constant segments \(y = \pm 1\):
This gives a single continuous nondecreasing solution over the entire real line.
(b) \(y(0)=1\): The constant function \(y(x)\equiv1\) is an obvious solution. Alternatively, using \(\sin(x+\alpha)=1\) at \(x=0\) gives \(\alpha=\dfrac{\pi}{2}+2k\pi\). Choosing the increasing branch near \(x=0\),
To extend for all \(x\):
There exist infinitely many such piecewise continuous solutions, obtained by shifting these arcs horizontally.
Step 5. Interpretation and summary.
The general solution \(y = \sin(x+\alpha)\) applies where \(|y|<1\) and \(y' = \sqrt{1-y^2} \ge 0\).
The singular solutions \(y = \pm1\) correspond to the flat portions of the slope field where \(y'=0\).
The complete solution set consists of rising sine arcs smoothly connected to horizontal segments at \(y=\pm1\).
All solutions.
y(x) = (x+) on intervals where (x+)0, together with y1.
8.3.2Linear First-Order Equations
A linear first-order equation has the form
Integrating Factor Method:
Compute the integrating factor:
\begin{equation} \mu(x) = e^{\int P(x)\,dx}. \end{equation}Multiply through by \(\mu(x)\):
\begin{equation} \frac{d}{dx}\!\left(\mu(x)y\right) = \mu(x)Q(x). \end{equation}Integrate:
\begin{equation} \mu(x)y = \int \mu(x)Q(x)\,dx. \end{equation}Solve for \(y(x)\):
\begin{equation} \boxed{y = [\mu(x)]^{-1} \int \mu(x) Q(x)\,dx}. \end{equation}
Solution.
Step 1. Identify \(P(x)\) and \(Q(x)\). The equation is in the standard linear form
where
Step 2. Compute the integrating factor. The integrating factor is
Step 3. Apply the general solution formula.
Substitute \(\mu(x) = e^{2x}\) and \(Q(x) = e^{-x}\):
Step 4. Evaluate the integral.
Final Answer:
Verification. Substitute \(y = e^{-x} + C e^{-2x}\) into the original equation:
which satisfies the ODE.
Separable equations and linear equations are the two types of first-order equations most commonly encountered in basic applications. However, there exist several other important forms of first-order equations that can be transformed into solvable types. In this section, we briefly introduce two such forms: the Bernoulli equation and exact equations. More details and applications are found in advanced differential equations texts.
8.3.4The Bernoulli Equation
The differential equation
where \(P(x)\) and \(Q(x)\) are known functions of \(x\), is called the Bernoulli equation. It is nonlinear for \(n \neq 0, 1\), but it can be transformed into a linear equation through an appropriate substitution.
Step 1. Multiply by \(y^{-n}\).
Step 2. Simplify using the substitution \(z = y^{1-n}\).
Hence, \(y^{-n}y' = \dfrac{1}{1-n}\,z'\). Substitute this into the previous equation:
Step 3. Multiply through by \((1-n)\) to simplify.
This is now a first-order linear ODE in \(z(x)\), which can be solved using the integrating factor method.
Reduced linear form: z' + (1-n)P(x)z = (1-n)Q(x).
After solving for \(z(x)\), the original dependent variable is recovered via
The Bernoulli equation provides an important link between nonlinear and linear first-order ODEs. By multiplying through by \(y^{-n}\) and substituting \(z = y^{1-n}\), we convert a nonlinear equation into a solvable linear form. This transformation frequently appears in problems involving population dynamics, exponential growth and decay, and fluid motion.
Solution.
Step 1. Multiply by \(y^{-n}\) with \(n=\tfrac{2}{3}\).
Step 2. Let \(z = y^{1-n}=y^{1/3}\). Then \(z'=\tfrac{1}{3}y^{-2/3}y'\), so \(y^{-2/3}y'=3z'\). Substitute:
Step 3. Solve the linear equation using the integrating-factor formula. For
the integrating factor is \(\mu(x)=e^{\int p(x)\,dx}=e^{x/3}\). Apply
Step 4. Back–substitute to \(y\). Since \(z=y^{1/3}\),
Step 5. Singular solution. Because \(y^{2/3}\ge 0\) for all real \(y\), the constant function
also satisfies the original ODE. Thus the complete set consists of the one-parameter family above together with \(y\equiv 0\).
8.3.5Exact Differential Equations
A general first-order differential equation can often be written in the form
where \(M\) and \(N\) are functions of two variables \(x\) and \(y\).
The first–order differential equation
is called exact if there exists a scalar function \(U(x,y)\) (called a potential function) such that
In this case, the general solution can be written implicitly as
A necessary and sufficient condition for exactness (assuming continuous partial derivatives) is
The condition above ensures that the mixed second derivatives of \(U(x,y)\) are equal:
This is guaranteed by Clairaut’s theorem when \(M\) and \(N\) have continuous partial derivatives.
Solution Procedure:
Check whether \( \dfrac{\partial M}{\partial y} = \dfrac{\partial N}{\partial x} \). If true, the equation is exact.
Integrate \(M(x,y)\) with respect to \(x\), treating \(y\) as constant:
\begin{equation} U(x,y) = \int M(x,y)\,dx + F(y), \end{equation}where \(F(y)\) is an arbitrary function of \(y\).
Since \(U(x,y)\) must satisfy
\begin{equation} \boxed{\,N(x,y) = \frac{\partial U}{\partial y}\,}, \end{equation}we differentiate the above expression for \(U(x,y)\) with respect to \(y\):
\begin{equation} \frac{\partial U}{\partial y} = \frac{\partial}{\partial y}\!\left(\int M(x,y)\,dx\right) + F'(y). \end{equation}Now substitute the known expression for \(N(x,y)\) from the original equation and equate:
\begin{equation} N(x,y) = \frac{\partial}{\partial y}\!\left(\int M(x,y)\,dx\right) + F'(y). \end{equation}From this, we can solve for \(F'(y)\):
\begin{equation} F'(y) = N(x,y) - \frac{\partial}{\partial y}\!\left(\int M(x,y)\,dx\right). \end{equation}Integrate with respect to \(y\) to find \(F(y)\), and then write the final expression for \(U(x,y)\). The general solution is obtained from
\begin{equation} \boxed{\,U(x,y) = \text{constant.}\,} \end{equation}Integrate \(F'(y)\) to find \(F(y)\), substitute back, and write the final solution \(U(x,y) = C\).
Solve
Solution.
Step 1. Identify
Step 2. Check the exactness condition:
They are equal, so the equation is exact.
Step 3. Integrate \(M(x,y)\) with respect to \(x\):
Step 4. Differentiate with respect to \(y\):
Compare with \(N(x,y) = x^2 + 3x y^2\), which implies \(F'(y)=0\).
Step 5. Hence \(F(y)=\text{constant}\). The potential function is
Step 6. The general solution is
8.3.6Integrating Factors
If the equation \(M(x,y)\,dx + N(x,y)\,dy = 0\) is not exact, it may sometimes be made exact by multiplying by a function \(\mu(x,y)\) called an integrating factor:
The new equation is exact if
Unfortunately, there is no universal method for finding \(\mu(x,y)\). However, there are two common and important special cases.
Case 1: Integrating factor depends on \(x\) only.
If \(\mu = \mu(x)\), then the equation becomes exact provided that
Hence,
Case 2: Integrating factor depends on \(y\) only.
If \(\mu = \mu(y)\), then the equation becomes exact provided that
Hence,
Solve
Solution.
Step 1. Identify \(M,N\) and test exactness.
Compute
Since \(M_y\neq N_x\), the equation is not exact.
Step 2. Seek an integrating factor \(\mu=\mu(x)\). Use
Hence
Step 3. Multiply the equation by \(\mu(x)=e^x\).
Define
Step 4. Check exactness \& find a potential \(U(x,y)\).
Integrate \(\tilde M\) w.r.t. \(x\) (treating \(y\) constant):
Using \(\int e^x(x+1)\,dx = x e^x\) and \(\int x e^x\,dx=(x-1)e^x\),
Now
Step 5. Implicit solution.
Solve
Solution.
Step 1. Identify \(M,N\) and test exactness.
Compute
Since \(M_y\neq N_x\), the equation is not exact.
Step 2. Seek an integrating factor \(\mu=\mu(y)\). Use
Hence
Step 3. Multiply the equation by \(\mu(y)=y^2\).
Define
Step 4. Check exactness \& find a potential \(U(x,y)\).
Integrate \(\tilde M\) w.r.t. \(x\) (treating \(y\) constant):
Now
Step 5. Implicit solution.
Many physical systems—such as thermodynamics, potential flow, and electrostatics—lead to exact differential equations, where \(U(x,y)\) represents a potential energy, enthalpy, or potential function. When the equation is not exact, an integrating factor often restores this structure.
8.4Second-Order Ordinary Differential Equations
A second-order linear ordinary differential equation (ODE) has the general form
where \(a_2(x), a_1(x), a_0(x)\), and \(f(x)\) are known functions of \(x\), and \(y(x)\) is the unknown function to be determined.
When \(a_2,a_1,a_0\) are constants, the equation simplifies to
We will first study the homogeneous case (\(f(x)=0\)) and then the nonhomogeneous case (\(f(x)\neq0\)).
8.4.1Homogeneous Equations with Constant Coefficients
A second–order homogeneous linear differential equation with constant coefficients has the form
Consider the homogeneous linear differential equation
Its solutions are determined by the characteristic (auxiliary) equation
If \(r_1\) and \(r_2\) are the roots of this quadratic, then the general solution is:
where \(C_1, C_2\) are arbitrary constants.
Consider the homogeneous constant–coefficient ODE
Step 1: Exponential trial and substitution. Assume a trial solution \(y=e^{r x}\) with constant \(r\). Then
Substituting into the ODE gives
Since \(e^{r x}\neq 0\) for all \(x\), we must have the characteristic equation
Step 2: Distinct real roots \(r_1\neq r_2\). If the characteristic polynomial has two distinct real roots \(r_1\) and \(r_2\), then
are both solutions. Their Wronskian is
so \(y_1,y_2\) are linearly independent. Hence the general solution is
Step 3: Repeated real root \(r_1=r_2=r\). (Construction of the \(x e^{r x}\) term.) If the characteristic polynomial has a double root \(r\), then we already have one solution \(y_1=e^{r x}\). To find a second, linearly independent solution, set
where \(v\) is to be determined. Compute derivatives:
Substitute \(y_2,y_2',y_2''\) into the ODE:
Divide by \(e^{r x}\neq 0\) and group terms by \(v'', v', v\):
Because \(r\) is a root of the characteristic polynomial, we have
Moreover, because the root is repeated, the derivative of the characteristic polynomial also vanishes at \(r\):
Hence the equation for \(v\) reduces to
Discarding the multiple of the first solution (the \(B e^{r x}\) part), we take \(v(x)=x\) and obtain
Therefore, the general solution in the repeated-root case is
Step 4: Complex conjugate roots \(r=\alpha\pm i\beta\) with \(\beta\neq 0\). If \(r=\alpha\pm i\beta\), then \(e^{(\alpha\pm i\beta)x}\) are solutions. Taking real and imaginary parts,
are real solutions. Their Wronskian is
so they are linearly independent. Thus the general real solution is
Combining the three cases completes the proof.
The step \(2ar+b=0\) in the repeated-root case is precisely the condition \(p'(r)=0\) where \(p(r)=ar^2+br+c\). This is why the reduction–of–order equation collapses to \(v''=0\), giving \(v(x)=x\) and hence \(y_2=x e^{r x}\).
The Wronskians shown above explicitly verify linear independence in the distinct-root and complex-root cases.
The exponential ansatz works here because constant coefficients make the ODE invariant under differentiation of \(e^{r x}\): derivatives only pull down factors of \(r\), converting the differential equation to an algebraic one.
8.4.2Nonhomogeneous Equations
A nonhomogeneous second-order linear ODE is
where \(f(x)\) is a known forcing (nonzero) function.
The general solution of a nonhomogeneous linear ODE is
where
\(y_h(x)\) is the homogeneous (complementary) solution satisfying \(a y''+b y'+c y=0\),
\(y_p(x)\) is a particular solution satisfying the full equation for the given \(f(x)\).
If \(f(x)\) is composed of elementary functions such as exponentials, polynomials, or sines and cosines, we assume a trial form for \(y_p(x)\) resembling \(f(x)\), substitute into the ODE, and solve for the unknown coefficients.
If any term of the assumed \(y_p(x)\) duplicates a part of \(y_h(x)\), multiply by the smallest power of \(x\) necessary to ensure linear independence.
8.4.3Example: Forced and Damped Harmonic Oscillator
Consider the driven, damped oscillator
where \(m\) is the mass, \(b\) the damping coefficient, \(k\) the spring constant, and \(F_0\sin(\omega t)\) the external driving force.
Solution.
Divide through by \(m\):
where
Step 2. Homogeneous solution.
For the unforced system,
the characteristic equation is
where
is the damped natural frequency. Hence,
Step 3. Particular (steady–state) solution using complex amplitude.
We look for a steady–state oscillation at the same frequency as the driving force. To simplify the algebra, we use a complex trial solution.
Step 3(a). Choose a complex trial solution.
Assume
where \(z_0\) is a complex constant. We may write
so that
Step 3(b). Substitute into the differential equation.
The driven damped oscillator equation (after dividing by \(m\)) is
We substitute the full complex function \(z_0 e^{i\omega t}\):
Substituting into the ODE gives
Step 3(c). Cancel \(e^{i\omega t}\).
Step 3(d). Solve for the complex amplitude \(z_0\).
Write the denominator in polar form:
where
Thus,
Step 3(e). Take the imaginary part.
Final steady–state response:
with
Step 4. Complete solution.
Step 5. Physical interpretation.
\(x_h(t)\): Transient response, which decays exponentially as \(e^{-\gamma t}\) and oscillates with the damped frequency \(\omega_d\).
\(x_p(t)\): Steady–state response, which oscillates at the driving frequency \(\omega\).
For weak damping (\(\gamma \ll \omega_0\)), the steady–state amplitude is
\begin{equation} A(\omega) = \frac{F_0/m}{\sqrt{(\omega_0^{\,2}-\omega^{2})^{2} + (2\gamma\omega)^{2}}}. \end{equation}The amplitude reaches a maximum (resonance) when the driving frequency is close to the natural frequency of oscillation,
\begin{equation} \omega \approx \omega_0. \end{equation}More precisely, minimising the denominator \((\omega_0^2-\omega^2)^2+(2\gamma\omega)^2\) gives the peak at \(\omega_{\text{res}} = \sqrt{\omega_0^2 - 2\gamma^2}\), which reduces to \(\omega_0\) in the weak-damping limit \(\gamma \ll \omega_0\). Evaluating the amplitude at \(\omega=\omega_0\),
\begin{equation} A(\omega_0) = \frac{F_0/m}{2\gamma\omega_0} = \frac{F_0}{2m\gamma\omega_0}, \end{equation}which is large when \(\gamma\) is small — the amplitude at resonance is limited only by the damping, and would diverge if \(\gamma\to0\).
As \(t \to \infty\), the transient term \(x_h(t)\) vanishes and only the steady oscillation \(x_p(t)\) remains.
8.5Where Differential Equations Are Used
A differential equation is a statement about how a quantity changes, and change is what physics is about. The forced damped oscillator worked at the end of this chapter is not one example among many — it is the template. The same equation \(m\ddot x + b\dot x + kx = F(t)\), with the symbols reinterpreted, governs a startling range of systems.
Mechanical vibration. Mass on a spring, a pendulum, a building swaying in wind: all are \(m\ddot x + b\dot x + kx = F(t)\), and the resonance condition \(\omega\approx\omega_0\) is why soldiers break step on bridges.
Electrical circuits. The series RLC circuit obeys \(L\ddot q + R\dot q + q/C = V(t)\), identical in form to the oscillator with \(L\leftrightarrow m\), \(R\leftrightarrow b\), \(1/C\leftrightarrow k\). Resonance here is how a radio selects a station.
Growth and decay. The first-order equation \(\dot y = ky\) describes radioactive decay, population growth, cooling, and the charging of a capacitor — the separable equations of this chapter.
Chemical kinetics and mixing. Reaction rates and concentration in a stirred tank are first-order linear ODEs solved by the integrating factor.
The starting point for everything later. Series solutions (Chapter 9) handle ODEs with variable coefficients that these elementary methods cannot; partial differential equations (Chapter 10) are solved by separating them into ordinary ones.
Two ideas from this chapter carry the most weight later. First, for a linear equation the general solution is \(y_h + y_p\): the full space of solutions is a particular solution plus the entire solution set of the homogeneous equation. Second, the exponential ansatz \(y=e^{rx}\) turns a constant-coefficient linear ODE into an algebraic equation, because differentiating \(e^{rx}\) only multiplies it by \(r\). Almost every linear system in physics is attacked by some version of these two moves.
Summary Table
| Type | Form | Method |
| Separable | \(\dfrac{dy}{dx} = g(x)h(y)\) | Separate and integrate: \(\int \dfrac{dy}{h} = \int g\,dx\). |
| Linear (1st order) | \(y' + P(x)y = Q(x)\) | Integrating factor \(\mu = e^{\int P\,dx}\). |
| Bernoulli | \(y' + P y = Q\,y^n\) | Substitute \(z = y^{1-n}\) to linearise. |
| Exact | \(M\,dx + N\,dy = 0\), \(M_y = N_x\) | Find \(U\) with \(U_x=M\), \(U_y=N\); solution \(U=C\). |
| Non-exact | \(M\,dx + N\,dy = 0\), \(M_y \neq N_x\) | Integrating factor \(\mu(x)\) or \(\mu(y)\). |
| Homogeneous (2nd, const.) | \(ay''+by'+cy = 0\) | Roots of \(ar^2+br+c=0\); three cases. |
| distinct real \(r_1,r_2\) | \(y = C_1 e^{r_1 x} + C_2 e^{r_2 x}\). | |
| repeated root \(r\) | \(y = (C_1 + C_2 x)e^{rx}\). | |
| complex \(\alpha\pm i\beta\) | \(y = e^{\alpha x}(C_1\cos\beta x + C_2\sin\beta x)\). | |
| Nonhomogeneous (2nd) | \(ay''+by'+cy = f(x)\) | \(y = y_h + y_p\); undetermined coefficients for \(y_p\). |
For the Interested Reader
Differential equations is a subject where the mechanics of solving can hide the point, which is that an equation relating a quantity to its rate of change already contains the system's whole future. The first item below is the best short antidote to losing sight of that. Everything listed is free.
The idea, seen
3Blue1Brown's opening film treats differential equations as a way of thinking about change, using the pendulum and the very oscillator of this chapter. It is worth watching before the techniques start to feel mechanical:
Differential equations, a tourist's guide — 3Blue1Brown
Watch on YouTube
Worked technique and practice
Paul's Online Math Notes — Differential Equations. The closest match to this chapter and the best source of extra practice: the full set, with Linear Equations and Second Order Equations covering our first- and second-order material in more detail, including variation of parameters, which we did not treat.
OpenStax, Calculus Volume 2, Chapter 4. 4.3 Separable Equations and 4.5 First-Order Linear Equations — a slower first pass with many worked examples and applications to growth, decay and mixing.
One method this chapter did not cover is worth knowing exists: the Laplace transform, which turns a linear constant-coefficient ODE with given initial conditions into an algebraic equation, handles discontinuous and impulsive forcing cleanly, and is the standard tool in engineering. Paul's notes above develop it in full. It is, in spirit, the same trade we made with \(e^{rx}\) and with the complex amplitude in the oscillator example — replace calculus by algebra by moving to a transform.