Title: Differential Equations

Q: Derivative
A: The \emph{derivative} of a function $f(x)$ with respect to $x$, interpreted as the rate of change of $f(x)$ with $x$, is
\begin{align*}
\frac{\mathrm{d} f}{\mathrm{d} x} = \lim_{h\to 0} \frac{f(x + h) - f(x)}{h}
\end{align*}
Alternatively:
\begin{align*}
\left.\frac{\mathrm{d} f}{\mathrm{d} x}\right\vert _{x_0} = \frac{f(x_0 + h) - f(x_0)}{h}+\frac{o(h)}{h}
\end{align*}
\begin{proof}
Let
\begin{align*}
\left.\frac{\mathrm{d} f}{\mathrm{d} x}\right\vert _{x_0} = \frac{f(x_0 + h) - f(x_0)}{h}+\frac{R(h)}{h}
\end{align*}
for some remainder functio $R(h)$. Now taking the limit of both sides as $h\to 0$, the LHS is a constant and the RHS has a term which can be recognised as the derivative:
\begin{align*}
\lim_{h\to 0}\left.\frac{\mathrm{d} f}{\mathrm{d} x}\right\vert _{x_0} &=\lim_{h\to 0}\left( \frac{f(x_0 + h) - f(x_0)}{h}+\frac{R(h)}{h}\right) \\
&=\lim_{h\to 0}\frac{f(x_0 + h) - f(x_0)}{h}+\lim_{h\to 0}\frac{R(h)}{h} \\
\therefore \left.\frac{\mathrm{d} f}{\mathrm{d} x}\right\vert _{x_0}&=\left.\frac{\mathrm{d} f}{\mathrm{d} x}\right\vert _{x_0} +\lim_{h\to 0}\frac{R(h)}{h}\\
\therefore \lim_{h\to 0}\frac{R(h)}{h}&=0 \\
\therefore R(h)&=o(h)
\end{align*}
Conversely:
\begin{align*}
\left.\frac{\mathrm{d} f}{\mathrm{d} x}\right\vert _{x_0} &= \frac{f(x_0 + h) - f(x_0)}{h}+\frac{o(h)}{h} \\
\left.\implies \frac{\mathrm{d} f}{\mathrm{d} x}\right|_{x_0} &= \lim_{h\to 0} \frac{f(x_0 + h) - f(x_0)}{h}
\end{align*}
\end{proof}

Q: Differentiable
A: For the function $f(x)$ to be differentiable, and so for the function $\frac{\mathrm{d} f}{\mathrm{d} x}$
to be well-defined at the point $x_0$, the left-hand limit (i.e. h is negative
and approaches zero from below) and the right-hand limit (h is positive and
approaches zero from above) must be defined and equal, i.e.
\begin{align*}
\lim_{h\to 0^-} \frac{f(x_0 + h) - f(x_0)}{h}=\lim_{h\to 0^+} \frac{f(x_0 + h) - f(x_0)}{h}
\end{align*}

Q: Order Parameters: Definitions
A: $f(x)=o[g(x)]$ as $x\to x_0$ if
\begin{align*}
\lim_{x\to x_0} \frac{f(x)}{g(x)}=0
\end{align*}
and \\
$F(x)=O[G(x)]$ as $x\to x_0$ if for all $x$ with $|x-x_0|<\delta$ (or if $x_0=\infty$, then for all $x>X$):
\begin{align*}
|F(x)|\leqslant M|G(x)|
\end{align*}
i.e. $\frac{f(x)}{g(x)}$ is bounded as $x\to x_0$

Q: Relationship between order parameters
A: \begin{align*}
f(x)=o[g(x)] \Rightarrow f(x)=O[g(x)]
\end{align*}
but \\
\begin{align*}
f(x)=O[g(x)] \nRightarrow f(x)=o[g(x)]
\end{align*}
as this is only true in the case $M=0$.

Q: Relationship between derivative and tangent line
A: \begin{align*}
\left.\frac{\mathrm{d} f}{\mathrm{d} x}\right\vert _{x_0} = \frac{f(x_0 + h) - f(x_0)}{h}+\frac{o(h)}{h}
\end{align*}
Multiplying by $h$:
\begin{align*}
h\left.\frac{\mathrm{d} f}{\mathrm{d} x}\right\vert _{x_0} = f(x_0 + h) - f(x_0)+o(h) \\
\therefore f(x_0 + h)=h\left.\frac{\mathrm{d} f}{\mathrm{d} x}\right\vert _{x_0}+f(x_0)+o(h) \\
\therefore y=(x-x_0)m+y_0+ o(h)
\end{align*}

Q: Chain rule
A: Let $f(x)=F[g(x)]$. Thus, from the definition of the derivative with repeated use of the tangent approximation equation:
\begin{align*}
\frac{\mathrm{d} f}{\mathrm{d} x} &= \lim_{h\to 0}\frac{f(x+ h) - f(x)}{h} \\
&= \lim_{h\to 0}\frac{F[g(x+ h)] - F[g(x)]}{h} \\
&= \lim_{h\to 0}\frac{F[g+h\frac{\mathrm{d} g}{\mathrm{d} x}+o(h)] - F[g]}{h} \\
&= \lim_{h\to 0}\frac{F(g+H) - F(g)}{h} \\
&= \lim_{h\to 0}\frac{\left(F(g)+H\frac{\mathrm{d} F}{\mathrm{d} g}+o(H)\right) - F(g)}{h} \\
&= \lim_{h\to 0}\frac{H\frac{\mathrm{d} F}{\mathrm{d} g}+o(H)}{h} \\
&= \lim_{h\to 0}\frac{\left(h\frac{\mathrm{d} g}{\mathrm{d} x}+o(h)\right)\frac{\mathrm{d} F}{\mathrm{d} g}+o(H)}{h} \\
&= \lim_{h\to 0}\frac{h\frac{\mathrm{d} g}{\mathrm{d} x}\frac{\mathrm{d} F}{\mathrm{d} g}+o(h)+o(H)}{h} \\
&= \frac{\mathrm{d} g}{\mathrm{d} x}\frac{\mathrm{d} F}{\mathrm{d} g}+\lim_{h\to 0}\frac{o(H)}{h}
\end{align*}
where $H=h\frac{\mathrm{d} g}{\mathrm{d} x}+o(h)$. We just need to check that $o(H)$ is $o(h)$ so the right hand term is zero. If $\frac{\mathrm{d} g}{\mathrm{d} x}=0$, then $H=o(h)$, so $o(H)$ is definitely $o(h)$. If $\frac{\mathrm{d} g}{\mathrm{d} x}\neq 0$, then we can still say $H=O(h)$, so $o(H)$ is $o(h)$.  \\
\\
Therefore, making the final observation that $F(g)$ is the same function as $f(g)$:
\begin{align*}
\frac{\mathrm{d} F[g(x)]}{\mathrm{d} x}=\frac{\mathrm{d} F}{\mathrm{d} g}\frac{\mathrm{d} g}{\mathrm{d} x}=\frac{\mathrm{d} f}{\mathrm{d} g}\frac{\mathrm{d} g}{\mathrm{d} x} \qed
\end{align*}

Q: Product Rule
A: Let $f(x)=u(x)v(x)$. Thus
\begin{align*}
\frac{\mathrm{d} f}{\mathrm{d} x} &= \lim_{h\to 0}\frac{f(x+ h) - f(x)}{h} \\
&= \lim_{h\to 0}\frac{u(x+ h)v(x+h) - u(x)v(x)}{h} \\
&= \lim_{h\to 0}\frac{\left(u+h\frac{\mathrm{d} u}{\mathrm{d} x}+o(h)\right)\left(v+h\frac{\mathrm{d} v}{\mathrm{d} x}+o(h)\right) - uv}{h} \\
&= \lim_{h\to 0}\frac{h\left( u\frac{\mathrm{d} v}{\mathrm{d} x}+v\frac{\mathrm{d} u}{\mathrm{d} x}\right) +h^2\frac{\mathrm{d} u}{\mathrm{d} x}\frac{\mathrm{d} v}{\mathrm{d} x}+o(h)}{h} \\
&= u\frac{\mathrm{d} v}{\mathrm{d} x}+v\frac{\mathrm{d} u}{\mathrm{d} x}+\lim_{h\to 0}h\frac{\mathrm{d} u}{\mathrm{d} x}\frac{\mathrm{d} v}{\mathrm{d} x}+\lim_{h\to 0}\frac{o(h)}{h} \\
&= u\frac{\mathrm{d} v}{\mathrm{d} x}+v\frac{\mathrm{d} u}{\mathrm{d} x} \qed
\end{align*}

Q: Quotient Rule
A: Letting $v\to v^{-1}$ in the product rule:
\begin{align*}
\frac{\mathrm{d}}{\mathrm{d} x}\left(\frac{u}{v}\right)&=v^{-1}\frac{\mathrm{d}u}{\mathrm{d} x}+u\left( -\frac{1}{v^2}\frac{\mathrm{d}v}{\mathrm{d} x}\right)\\
&=\frac{v\frac{\mathrm{d}u}{\mathrm{d} x}-u\frac{\mathrm{d}v}{\mathrm{d} x}}{v^2} \qed
\end{align*}
Alternatively, let's go back to first principles. Let $f(x)=\frac{u(x)}{v(x)}$. Thus
\begin{align*}
\frac{\mathrm{d} f}{\mathrm{d} x} &= \lim_{h\to 0}\frac{f(x+ h) - f(x)}{h} \\
&= \lim_{h\to 0}\frac{\frac{u(x+ h)}{v(x+h)} - \frac{u(x)}{v(x)}}{h} \\
&= \lim_{h\to 0}\frac{u(x+ h)v(x) - u(x)v(x+h)}{hv(x+h)v(x)} \\
&= \lim_{h\to 0}\frac{v(u+h\frac{\mathrm{d}u}{\mathrm{d} x}+o(h))-u(v+h\frac{\mathrm{d}v}{\mathrm{d} x}+o(h))}{hv(v+h\frac{\mathrm{d}v}{\mathrm{d} x}+o(h))} \\
&= \lim_{h\to 0}\frac{h(v\frac{\mathrm{d}u}{\mathrm{d} x}-u\frac{\mathrm{d}v}{\mathrm{d} x})+o(h)}{hv(v+h\frac{\mathrm{d}v}{\mathrm{d} x}+o(h))} \\
&= \lim_{h\to 0}\frac{v\frac{\mathrm{d}u}{\mathrm{d} x}-u\frac{\mathrm{d}v}{\mathrm{d} x}+\frac{o(h)}{h}}{v(v+h\frac{\mathrm{d}v}{\mathrm{d} x})+o(h)}\\
&=\frac{v\frac{\mathrm{d}u}{\mathrm{d} x}-u\frac{\mathrm{d}v}{\mathrm{d} x}}{v^2} \qed
\end{align*}

Q: (General) Leibniz rule
A: \begin{align*}
[uv]^{(n)}=\sum_{r=0}^n {n\choose r} u^{(r)}v^{(n-r)}
\end{align*}
\begin{proof}
Let $\varphi(n)$ be the statement
\begin{align*}
[uv]^{(n)}=\sum_{r=0}^n {n\choose r} u^{(r)}v^{(n-r)}
\end{align*}
Clearly $\varphi(0)$ is trivial and $\varphi(1)$ is just the product rule. \\
\\
Suppose $\exists k\geqslant 0: \varphi(k)$.
\begin{align*}
[uv]^{(k)}=\sum_{r=0}^k {k\choose r} u^{(r)}v^{(k-r)}
\end{align*}
Thus, differentiating and considering the coefficient of $u^{(r)}$:
\begin{align*}
[uv]^{(k+1)}&={k\choose 0}u^{(0)}v^{(k+1)}+{k+1\choose 0}u^{(k+1)}v^{(0)} +\sum_{r=1}^{k} u^{(r)}\left(v^{(k-r+1)}{k\choose r}+ v^{(k-(r-1))}{k\choose r-1}\right) \\
&={k\choose 0}u^{(0)}v^{(k+1)}+{k+1\choose 0}u^{(k+1)}v^{(0)} +\sum_{r=1}^{k} u^{(r)}v^{(k+1-r)}\left({k\choose r}+ {k\choose r-1}\right) \\
&={k+1\choose 0}u^{(0)}v^{(k+1)}+{k+1\choose 0}u^{(k+1)}v^{(0)} +\sum_{r=1}^{k} {k+1\choose r}u^{(r)}v^{(k+1-r)}\\
&=\sum_{r=0}^{k+1} {k+1\choose r}u^{(r)}v^{(k+1-r)}
\end{align*}
Thus $\varphi(k)\implies \varphi(k+1)$ for any $k\geqslant 1$. So by induction, $\forall n\geqslant 0: \varphi(n)$.
\end{proof}

Q: Taylor's Theorem
A: Suppose
\begin{align*}
f(x)=\sum_{r=0}^\infty a_r(x-x_0)^r.
\end{align*}
Then, $f^{(r)}(x_0)=r!a_r$. Thus
\begin{align*}
f(x)=f(x_0)+(x-x_0)\left.\frac{\mathrm{d}f}{\mathrm{d}x}\right\vert_{x=x_0}+\frac{(x-x_0)^2}{2!}\left.\frac{\mathrm{d}^2f}{\mathrm{d}x^2}\right\vert_{x=x_0}\ldots + \frac{(x-x_0)^n}{n!}\left.\frac{\mathrm{d}^nf}{\mathrm{d}x^n}\right\vert_{x=x_0} +o((x-x_0)^n)
\end{align*}
or
\begin{align*}
f(x_0+h)=f(x_0)+h\left.\frac{\mathrm{d}f}{\mathrm{d}x}\right\vert_{x=x_0}+\frac{h^2}{2!}\left.\frac{\mathrm{d}^2f}{\mathrm{d}x^2}\right\vert_{x=x_0}\ldots + \frac{h^n}{n!}\left.\frac{\mathrm{d}^nf}{\mathrm{d}x^n}\right\vert_{x=x_0} +o(h^n)
\end{align*}

Q: Indeterminate forms
A: When replacing factors of a function by their limit does not give enough information to determine the limit of the function e.g. if $\frac{0}{0}$ or $\frac{\infty}{\infty}$ is obtained.

Q: L'H\^opital's Rule
A: \begin{align*}
\lim_{x\to x_0} \frac{f(x)}{g(x)}=\lim_{x\to x_0} \frac{f'(x)}{g'(x)}
\end{align*}
provided $g'(x_0)\neq 0$ and $f,g$ are differentiable at $x=x_0$. If $\frac{f'(x)}{g'(x)}$ is still an indeterminate form, then you can go to a higher derivative (by induction). \\
\\
Consider the special case $f(x_0)=0=g(x_0)$. Thus
\begin{align*}
\lim_{x\to x_0} \frac{f(x)}{g(x)}&=\lim_{x\to x_0} \frac{0+hf'+o(h)}{0+hg'+o(h)}\\
&=\lim_{x\to x_0} \frac{f'+\frac{o(h)}{h}}{g'+\frac{o(h)}{h}}\\
=\lim_{x\to x_0} \frac{f'}{g'} \qed
\end{align*}

Q: Riemann Sum
A: \begin{align*}
\int_a^b f(x) \mathrm{d} x = \lim_{\Delta x\to 0}\sum_{n=0}^{N-1} f(x_n)\Delta x=\lim_{N\to\infty}\sum_{n=0}^{N-1} f(x_n)\Delta x
\end{align*}
where
\begin{align*}
\Delta x =\frac{b-a}{N},x_n=a+n\Delta x
\end{align*}
Consider the rectangle with base $x_n\leqslant x\leqslant x_{n+1}$. By the MVT, for all $x$ there is a point $x_c$ such that $x_n<x_c<x$ and
\begin{align*}
f(x)-f(x_n)=(x-x_n)f'(x_c)
\end{align*}
which gives us the following upper bound on the height of the curve above the rectangle (in the Riemann sum approximation to the area under the curve):
\begin{align*}
\left| f(x)-f(x_n) \right| &\leqslant \left|(x_{n+1}-x_n)\right| \left|\max_{x\in[a,b]}f'(x_c)\right| \\
&=M\Delta x
\end{align*}
as we are considering the interval $x_n\leqslant x\leqslant x_{n+1}$ and we can find the largest $f'(x_c)$ for well-behaved functions (N.B. $M$ is a fixed positive number). \\
\\
Thus the error $\epsilon$ in the \textbf{area} is bounded above by $|\epsilon|\leqslant M(\Delta x)^2$ so $|\epsilon|=O((\Delta x)^2)=O\left(\frac{b-a}{N}\right)^2$. \\
\\
Summing these errors over the $N$ strips gives that the error in the \textbf{total} area is $O\left(\frac{(b-a)^2}{N}\right)$ which indeed goes to zero as $N\to\infty$ as required. \qed

Q: Mean Value Theorem
A: For well-behaved functions, for any interval $[a,b]$ there is a point $c \in [a,b]$:
\begin{align*}
f(b)=f(a)+(b-a)f'(c)
\end{align*}
i.e. there is always a point with gradient equal to that of the line AB. \\
\\
This is a generalised version of Rolle's Theorem which assumes that $f(a)=f(b)$ (the two are equivalent by turning the page).

Q: Fundamental Theorem of Calculus
A: \begin{align*}
F(x)=F(a)+\int_a^xf(t)\mathrm{d}t \iff \frac{\mathrm{d}F}{\mathrm{d}x}=f(x)
\end{align*}
\begin{proof}
By the definition of a derivative
\begin{align*}
\frac{\mathrm{d}F}{\mathrm{d}x}&=\lim_{h\to 0}\frac{1}{h}  \left[ \int_a^{x+h}f(t)\mathrm{d}t-\int_a^xf(t)\mathrm{d}t \right] \\
&=\lim_{h\to 0}\frac{1}{h} \int_x^{x+h}f(t)\mathrm{d}t \\
&=\lim_{h\to 0}\frac{1}{h}\left[ f(x)h+O(h^2)\right] \\
&=f(x)+\lim_{h\to 0}O(h) \\
&=f(x)
\end{align*}
The reverse implication follows by uniqueness of the solution of first order linear ODEs.
\end{proof}

Q: Differentiation under the integral sign
A: Define $I(b, c) = \int_0^bf(t, c)\mathrm{d} t$. Then
\begin{align*}
  \frac{\partial I}{\partial c} &= \lim_{\delta c\to 0}\frac{1}{\delta c}\left[\int_0^bf(t, c + \delta c)\mathrm{d} t - \int_0^b f(t, c)\mathrm{d} t\right]\\
  &= \lim_{\delta c\to 0}\int_0^b\frac{f(t, c + \delta c) - f(t, c)}{\delta c}\mathrm{d} t\\
  &= \int_0^b\lim_{\delta c\to 0}\frac{f(t, c + \delta c) - f(t, c)}{\delta c}\mathrm{d} t\\
  &= \int_0^b\frac{\partial f}{\partial c}\mathrm{d}t
\end{align*}
Now, using the multivariate chain rule:
\begin{equation*}
\frac{\mathrm{d}}{\mathrm{d}x}\int_{a(x)}^{b(x)} f(x,t) \ \mathrm{d}t=\\
f(x,b(x))b'(x)-f(x,a(x))a'(x)+\int_{a(x)}^{b(x)}\frac{\partial}{\partial x}f(x,t) \ \mathrm{d}t
\end{equation*}

Q: Integration by parts
A: From the product rule
\begin{align*}
\frac{\mathrm{d}}{\mathrm{d}x}uv&=uv'+vu' \\
\implies uv &= \int uv' \mathrm{d}x + \int vu' \mathrm{d}x \\
\implies \int uv' \mathrm{d}x &= uv- \int vu' \mathrm{d}x
\end{align*}

Q: Partial Derivative
A: Given a function of several variables $f(x, y)$, the \emph{partial derivative} of $f$ with respect to $x$ is the rate of change of $f$ as $x$ varies, keeping $y$ constant. It is given by
  \[
    \left. \frac{\partial f}{\partial x}\right|_y = \lim_{\delta x\to 0} \frac{f(x + \delta x, y) - f(x, y)}{\delta x}
  \]

Q: Schwarz's (Clairaut's) Theorem
A: Symmetry of second derivatives: $f_{xy}=f_{yx}$ where
\[
\frac{\partial ^2 f}{\partial x \partial y}=\frac{\partial }{\partial x}\frac{\partial f}{\partial y}=f_{yx}
\]
(holds if second derivatives are continuous)

Q: Multivariate Chain Rule
A: \begin{align*}
\delta f &= f(x+\delta x,y+\delta y)-f(x,y) \\
&=f(x+\delta x,y+\delta y)-f(x+\delta x,y)+f(x+\delta x,y)-f(x,y) \\
&= \frac{\partial f}{\partial y}(x+\delta x,y) \delta y +o(\delta y) + \frac{\partial f}{\partial x}(x,y) \delta x +o(\delta x) \\
&= \left[\frac{\partial f}{\partial y}(x,y)+ \delta x\frac{\partial^2 f}{\partial x \partial y}(x,y)+o(\delta x) \right]\delta y +o(\delta y) + \frac{\partial f}{\partial x}(x,y) \delta x +o(\delta x) \\
&= \frac{\partial f}{\partial y}\delta y+\frac{\partial f}{\partial x} \delta x + \frac{\partial^2 f}{\partial x \partial y}\delta x\delta y+o(\delta x) +o(\delta y) \\
&= \frac{\partial f}{\partial y}\delta y+\frac{\partial f}{\partial x} \delta x +o(\delta x) +o(\delta y)
\end{align*}
This gives the two forms of the multivariate chain rule when the limit is taken (known as differential and integral forms)
\begin{align*}
\frac{\mathrm{d} f}{\mathrm{d} t} &= \frac{\partial f}{\partial x}\frac{\mathrm{d} x}{\mathrm{d} t} + \frac{\partial f}{\partial y}\frac{\mathrm{d} y}{\mathrm{d} t} \\
\int \mathrm{d} f &= \int \frac{\partial f}{\partial x}\mathrm{d} x + \int \frac{\partial f}{\partial y}\mathrm{d} y
\end{align*}
This can be generalised to
\begin{align*}
\delta f = &f(x_1+\delta x,x_2+\delta x_2 \ldots x_{n-1}+\delta x_{n-1},x_n + \delta x_n)-f(x_1,x_2,\ldots x_{n-1},x_n) \\
=&f(x_1+\delta x,x_2+\delta x_2 \ldots x_{n-1}+\delta x_{n-1},x_n + \delta x_n)-f(x_1,x_2+\delta x_2 \ldots x_{n-1}+\delta x_{n-1}, x_n + \delta x_n) \\
+&f(x_1,x_2+\delta x_2 \ldots x_{n-1}+\delta x_{n-1}, x_n + \delta x_n)-f(x_1,x_2 \ldots x_{n-1}+\delta x_{n-1}, x_n + \delta x_n) \\
+&\ldots  \\
+&f(x_1,x_2 \ldots x_{n-1},x_n + \delta x_n)-f(x_1,x_2,\ldots x_{n-1},x_n) \\
=&\frac{\partial f}{\partial x_1}\mathrm{d}x_1+\frac{\partial f}{\partial x_2}\mathrm{d}x_2+\ldots +\frac{\partial f}{\partial x_{n-1}}\mathrm{d}x_{n-1}+\frac{\partial f}{\partial x_n}\mathrm{d}x_n
\end{align*}
N.B. You must take into account all the dimensions i.e. all the variables which could change $f$.

Q: Multivariate implicit differentiation
A: Given an equation
  \[
    F(x, y, z) = c
  \]
  for some constant $c$, we have
\begin{align*}
    \mathrm{d}  F &= \frac{\partial F}{\partial x}\mathrm{d} x + \frac{\partial F}{\partial y}\mathrm{d}  y + \frac{\partial F}{\partial z}\mathrm{d}  z\\
    \left.\frac{\partial F}{\partial x}\right|_y &= \frac{\partial F}{\partial x}\left.\frac{\partial x}{\partial x}\right|_y + \frac{\partial F}{\partial y}\left.\frac{\partial y}{\partial x}\right|_y + \frac{\partial F}{\partial z}\left.\frac{\partial z}{\partial x}\right|_y \\
   &=\frac{\partial F}{\partial x} + \frac{\partial F}{\partial z}\left.\frac{\partial z}{\partial x}\right|_y \\
&= 0\\
    \implies\left.\frac{\partial z}{\partial x}\right|_y &= -\frac{(\partial F)/(\partial x)}{(\partial F)/(\partial z)}\\
\therefore \left.\frac{\partial z}{\partial x}\right|_y \left.\frac{\partial x}{\partial y}\right|_z \left.\frac{\partial y}{\partial z}\right|_x &=-1
  \end{align*}

Q: The exponential function
A: Consider $f(x)=a^x$. The derivative of this function is given by
\begin{align*}
  \frac{\mathrm{d} f}{\mathrm{d} x} &= \lim_{h\to 0}\frac{a^{x+h}-a^x}{h}\\
  &= a^x\lim_{h\to 0}\frac{a^h-1}{h}\\
  &= \lambda a^x\\
  &= \lambda f(x)
\end{align*}
where $\displaystyle \lambda = \lim_{h\to 0}\frac{a^h-1}{h} = f'(0) = $ const. So the derivative of an exponential function is a multiple of the original function. This suggests there must be some $a$ for which the exponential function is its own derivative. \\
\\
This prompts the definition $\exp(x) = e^x$ as the unique function $f$ satisfying $f'(x) = f(x)$ and the initial condition $f(0) = 1$. It can be shown that \[
e=\lim_{k\to\infty} \left( 1+\frac{1}{k} \right)^k
\]
which gives the numerical value of $e$.

Q: Linear Differential Equation
A: The dependent variable only appears linearly in the equation

Q: Homogeneous Differential Equation
A: $y=0$ is a solution to the equation i.e. there is no forcing function in $x$

Q: Order of a differential equation
A: The highest derivative which appears in the equation

Q: Well-posed problem
A: Problem with initial or boundary conditions where you can find the specific solution.

Q: forward Euler approximation
A: This is a method of approximating a differential equation with a discrete system using the transformations
\begin{align*}
\frac{\mathrm{d}y}{\mathrm{d}x}=\frac{y_{n+1}-y_n}{h} \\
x=nh
\end{align*}

Q: Series solution of FOLHODE with CCs
A: Let
\begin{align*}
y&=\sum_{n=0}^\infty a_nx^n \\
\implies y'&=\sum_{n=1}^\infty na_nx^{n-1} \\
\end{align*}
so
\begin{align*}
by'+cy=0 \\
\implies bxy'+cxy=0 \\
\implies bx\sum_{n=1}^\infty na_nx^{n-1}+cx\sum_{n=0}^\infty a_nx^n=0 \\
\implies \sum_{n=1}^\infty \left( bna_n+ca_{n-1} \right)x^n =0 \\
\implies bna_n+ca_{n-1}=0 \\
\implies a_n=-\frac{ca_{n-1}}{bn} \\
\implies a_n=\left( -\frac{c}{b}\right)^n\frac{a_0}{n!} \\
\implies y=a_0 \sum_{n=0}^\infty  \frac{(-cx/b)^n}{n!} \\
\implies y=a_0 \exp\left(-\frac{cx}{b}\right)
\end{align*}
The last line shows the consistency with the eigenfunction and separation of variables solution

Q: Integrating factor
A: Consider the general form of equation
\[
  a(x)y' + b(x)y = c(x).
\]
Divide by $a(x)$ to get the standard form
\[
  y' + p(x) y = f(x).
\]
We solve this by multiplying an integrating factor $\mu (x)$ to obtain $(\mu)y' + (\mu p)y = \mu f$.
We want to choose a $\mu$ such that the left hand side is equal to $(\mu y)'$. By the product rule, we want $\mu p = \mu'$, ie.
\begin{align*}
  p &= \frac{1}{\mu}\frac{\mathrm{d} \mu }{\mathrm{d} x}\\
  \int p \mathrm{d} x &= \int \frac{1}{\mu}\frac{\mathrm{d}\mu}{\mathrm{d} x} \mathrm{d} x\\
  &= \int \frac{1}{\mu} \mathrm{d} u\\
  &= \ln \mu (+ C)\\
  \mu &= \exp\left(\int p \mathrm{d} x\right)
\end{align*}
Then by construction, we have $(\mu y)' = \mu f$ and thus
\[
  y = \frac{\int \mu f \mathrm{d} x}{\mu},\text{ where } \mu = \exp\left(\int p \mathrm{d} x\right)
\]

Q: Separable equations
A: A first-order differential equation is \emph{separable} if it can be manipulated into the following form:
  \[
    q(y) \mathrm{d} y = p(x) \mathrm{d} x.
  \]
  in which case the solution can be found by integration
  \[
    \int q(y)\mathrm{d} y = \int p(x)\mathrm{d} x
  \]

Q: Exact equation
A: The nonlinear homogeneous ODE
\[
  Q(x, y)\frac{\mathrm{d} y}{\mathrm{d} x} + P(x, y) = 0.
\]
is exact if there is a function $f$ such that
\[
\mathrm{d} f = Q(x, y)\mathrm{d} y + P(x, y)\mathrm{d} x
\]
 For simply connected domains, we can identify the condition for exactness as follows. Using the chain rule
\[
  \frac{\partial f}{\partial x} = P,\quad\frac{\partial f}{\partial y} = Q.
\]
From this we have
\[
  \frac{\partial^2 f}{\partial y\partial x} = \frac{\partial P}{\partial y},\quad\frac{\partial^2 f}{\partial x \partial y} = \frac{\partial Q}{\partial x}.
\]
We know that the two mixed 2nd derivatives are equal. So the condition for exactness in simply connected domains is
\[
  \frac{\partial P}{\partial y} = \frac{\partial Q}{\partial x}.
\]
If the equation is exact, then the solution is simply $f = $ constant, and we can find $f$ by integrating $\frac{\partial f}{\partial x} = P$ (giving you a $g(y)$ term), then differentiate the answer and compare to $\frac{\partial f}{\partial y} = Q$ (solving for $g(y)$ and leaving only a constant of integration).

Q: Simply-connected domain
A: A domain $\mathcal{D}$ is simply-connected if it is connected and any closed curve in $\mathcal{D}$ can be shrunk to a point in $\mathcal{D}$ without leaving $\mathcal{D}$.

Q: Isoclines
A: Given \[
\frac{\mathrm{d}y}{\mathrm{d}t}=f(y,t)
\]
we can construct curves of constant slope given by $f(y,t)=C$. These can then be connected appropriately by drawing a curve which has gradient $C$ crossing each isocline. \\
\\
Rule: for single-valued functions, the curves cannot cross.

Q: Fixed (equilibrium) points
A: A fixed point of \[
\frac{\mathrm{d}y}{\mathrm{d}t}=f(y,t)
\]
is a solution for which $\frac{\mathrm{d}y}{\mathrm{d}t}=0$ or equivalently $f(y,t)=0$ which gives a solution of the form $y=c$.

Q: Peturbation analysis
A: Consider the equation \[
\frac{\mathrm{d}y}{\mathrm{d}t}=f(y,t)
\]
Let $y_f$ be a fixed point. Now we consider a small perturbation $\epsilon (t)$ from this fixed point, such that $y(t)=y_f+\epsilon (t)$. Differentiating both sides and using the differential equation gives
\begin{align*}
\frac{\mathrm{d}y}{\mathrm{d}t}&=0+\frac{\mathrm{d}\epsilon }{\mathrm{d}t} \\
\implies\frac{\mathrm{d} \epsilon }{\mathrm{d}t}&=f(y,t) \\
&=f(y_f+\epsilon,t) \\
&=f(y_f,t)+\epsilon  \frac{\partial f}{\partial y} (y_f,t) +\frac{\epsilon ^2 }{2} \frac{\partial ^2 f}{\partial y^2} (y_f,t)+o(\epsilon ^2) \\
&=\epsilon  \frac{\partial f}{\partial y} (y_f,t) +\frac{\epsilon ^2 }{2} \frac{\partial ^2 f}{\partial y^2} (y_f,t)+o(\epsilon ^2) \\
\end{align*}
Thus the flow in the vicinity of the fixed point can be understood by the equation \[
\frac{\mathrm{d} \epsilon }{\mathrm{d}t}=\frac{\epsilon ^n }{n!} \frac{\partial ^n f}{\partial y^n} (y_f,t)
\]
where $n$ is the least value such that $\frac{\partial ^n f}{\partial y^n}(y_f,t)\neq 0$. You can then either solve explicitly for $\epsilon(t)$ and observe the behaviour as $t\to\infty$ or you can simply consider the sign of $\frac{\partial ^n f}{\partial y^n} (y_f,t)$ and the parity of $n$. Specifically, if $n$ is odd and $\frac{\partial ^n f}{\partial y^n} (y_f,t)$ is negative, then $y_f$ is stable, otherwise unstable.

Q: Autonomous system
A: A system governed by the equation \[
\frac{\mathrm{d}y}{\mathrm{d}t}=f(y)
\]
so has no explicit dependence on time. \\
\\
This system thus satisfies
\begin{align*}
\frac{\mathrm{d} \epsilon }{\mathrm{d}t}&=\epsilon  \frac{\partial f}{\partial y} (y_f,t) \\
&= \epsilon\frac{\mathrm{d} f}{\mathrm{d}y} (y_f) \\
&= k\epsilon\\
\end{align*}
so \[
\epsilon=\epsilon_0e^{kt}
\]
and thus $y_f$ is stable iff. $k<0$.

Q: Phase portrait
A: Draw the dependent variable on a horizontal axis. Where the gradient is positive, draw arrows pointing in the positive $x$ direction. Where the gradient is negative, draw arrows pointing in the negative $x$ direction. Where the gradient is zero, draw vertical lines. \\
\\
The arrows now point towards stable fixed points and away from unstable fixed points.

Q: Differential Logistic Equation
A: Now imagine that the environment is finite, and that there is a need for food. The birth rate is $\alpha$ and the death rate is $\beta$. The probability of some food being found by some individual member of the population may be modelled as being proportional to the population size $y$. Therefore, the probability that the same source of food is found by two
individuals is proportional to the square of the population size, i.e. $y^2$. If food is scarce,  then there will be a fight to the death at a rate $\gamma$.
\begin{align*}
  \frac{\mathrm{d} y}{\mathrm{d} t} &= (\alpha - \beta)y - \gamma y^2\\
  \frac{\mathrm{d} y}{\mathrm{d} t}&= \lambda y \left(1 - \frac{y}{Y}\right),
\end{align*}
where $\lambda = \alpha - \beta$ and $Y = \lambda / \gamma $. This is the differential logistic equation. Note that it is separable and can be solved explicitly. Note also that the population increases until it reaches its stable fixed point at $y=Y$.

Q: Logistic Map
A: The discrete version of the logistic equation is
\begin{align*}
  \frac{x_{n+1} - x_n}{\Delta t} &= \lambda x_n -\gamma x_n^2 \\
  x_{n + 1} &= (1+\lambda\Delta t)x_n-\gamma\Delta t x_n^2
\end{align*}
We assume the time interval is choosen so that $r=1+\lambda\Delta t=\gamma\Delta t$. Thus \[
 x_{n + 1}=rx_n(1-x_n)
\]

Q: Cobweb diagram
A: \begin{itemize}
\item Start from an initial value $x_0$ on the $x$-axis.
\item Draw a vertical line up to the plot of $f(x)$. This determines $x_1$.
\item Draw a horizontal line to the diagonal $y = x$. This is thus at the
appropriate $x = x_1$ value.
\item Draw a vertical line to the plot of $f(x)$ to determine $x_2 \ldots $
\item Repeat this sequence of: vertical line to $f(x)$ then horizontal line to
$y = x$ as often as required.
\end{itemize}
e.g. $f(x)=rx(1-x)$

Q: Period $k$ limit cycle
A: $x_{n+k}=x_n$ i.e. an oscillating iteration which does not converge to a fixed point

Q: Fixed points of logistic map
A: If $x_{n+1}=x_n$, then \[
x_n=rx_n(1-x_n)
\]
giving two solutions: $x_n=0$ and $x_n=1-\frac{1}{r}$. However, for $r < 1$, $1-\frac{1}{r}<0$ so we ignore this fixed point. \\
\\
Consider the fixed point $x_n=0$. Let the error of the $n$th iterate from the fixed point be $\epsilon_n$. Therefore \[
\epsilon_{n+1}=r\epsilon_n(1-\epsilon_n)\approx r\epsilon_n
\]
for sufficiently small $\epsilon_n$. Thus $x_n=0$ is a stable fixed point iff. $r<1$. \\
\\
Now consider $x_f=1-\frac{1}{r}$. Let $x_n=x_f+\epsilon_n$ as before. Therefore  \vspace{-0.5cm}
\begin{align*}
x_{n+1}&=rx_n(1-x_n) \\
x_f+\epsilon_{n+1}&=r (x_f+\epsilon_n)(1-[x_f+\epsilon_n]) \\
1-\frac{1}{r}+\epsilon_{n+1}&=r \left( 1-\frac{1}{r}+\epsilon_n\right) \left( 1-\left[1-\frac{1}{r}+\epsilon_n\right] \right) \\
&\approx 1-\frac{1}{r} +\epsilon_n-r\epsilon_n+\epsilon_n \\
\implies \epsilon_{n+1}&\approx (2-r)\epsilon_n
\end{align*}
Thus $x_f=1-\frac{1}{r}$ is stable iff. $1<r<3$ and in particular, $x_f$ exhibits oscillatory convergence iff. $2<r<3$. \\
\\
An alternate way to check the stability is using the Taylor series expansion of $f(x)=rx(1-x)$:
\begin{align*}
f(x_f+\epsilon )&=f(x_f)+\epsilon \frac{\mathrm{d}f}{\mathrm{d}x} (x_f)+o(\epsilon) \\
&=x_f+\epsilon \frac{\mathrm{d}f}{\mathrm{d}x} (x_f)+o(\epsilon ) \\
&\approx x_f+\epsilon \frac{\mathrm{d}f}{\mathrm{d}x} (x_f) \\
\end{align*}
Thus $x_f$ is stable iff.  \[
\left| \frac{\mathrm{d}f}{\mathrm{d}x} (x_f) \right| <1
\]
So in this case \[
\frac{\mathrm{d}f}{\mathrm{d}x}=r(1-2x)
\]
This gives the same pattern of signs as before for each fixed point.
\[
\left.\frac{\mathrm{d}f}{\mathrm{d}x}\right|_0=r
\]
or
\[
\left.\frac{\mathrm{d}f}{\mathrm{d}x}\right|_{1-1/r} =2-r
\]

Q: Bifurcations
A: As $r$ increases beyond $3$ in the logistic map equation, non-converging oscillations appear which can be understood as fixed points of the second iterate. As $r$ increases past $3.52$, a period $4$ limit cycle then appears (i.e. a fixed point of the fourth iterate). In fact, the distance between these bifurcations decreases with $r$ and the ratio between the successive distances converges to $\delta\approx 4.7$ (the Feigenbaum constant).

Q: Requirements for solutions to form a basis of the space of solutions of the homogeneous equation
A: The solutions must be linearly independent and complete to span the set

Q: Magical detuning technique for second order ODEs with constant coefficients
A: Given a 2nd order LHODE with a repeated root \[
y''-2ky'+k^2y=0
\]
with characteristic equation \[
(\lambda -k)^2=0
\]
we only get one solution, namely $y_1=e^{kx}$. In order to find a linearly independent solution, we consider the \emph{detuned} equation \[
y''-2ky'+(k^2-\epsilon^2)y=0
\]
with characteristic equation \[
(\lambda -k)^2 =\epsilon^2
\]
and roots $\lambda=k\pm\epsilon$. This gives the complementary function
\begin{align*}
y&=Ae^{(k+\epsilon)x}+Be^{(k-\epsilon)x} \\
&=e^{kx} [Ae^{\epsilon x}+Be^{-\epsilon x}] \\
&=e^{kx} [A(1+\epsilon x + O(\epsilon^2) )+B(1-\epsilon x + O(\epsilon^2))] \\
&=e^{kx} [(A+B)+\epsilon x(A-B)+O(A\epsilon^2)+O(B\epsilon^2)] \\
&=e^{kx} [\alpha+ \beta x+O(A\epsilon^2)+O(B\epsilon^2) ]
\end{align*}
where $\alpha=A+B$ and $\beta=\epsilon(A-B)$. This means that $A = \frac{1}{2}(\alpha + \frac{\beta}{\epsilon})$ and $B = \frac{1}{2}(\alpha - \frac{\beta}{\epsilon})$. We take $\alpha$ and $\beta$ as fixed constants independent of $\epsilon$, so $A$ and $B$ diverge as $\epsilon \to 0$ in. In fact, $A = O(\frac{1}{\epsilon})$, so we know that $O(A\epsilon^2) = O(\epsilon)$, and similarly for $B$. So
\begin{align*}
y&=e^{kx} [\alpha+ \beta x+O(\epsilon) ] \\
\implies\lim_{\epsilon\to 0} y(x,\epsilon)&=e^{kx} [\alpha+ \beta x]
\end{align*}
This very particular choice of limit was necessary to prevent $\beta$ from going to zero and getting the alternate, independent solution. \\
\\
This is a general rule for linear equations with constant coefficients. If $y_1(x)$ is a degenerate complementary function, then $y_2(x)=xy_1(x)$ is an independent complementary function.

Q: Reduction of order
A: Given the equation \[
y''+p(x)y'+q(x)y=0
\]
and a first complementary function $y_1(x)$, we can find an independent compementary function (in fact the entire complementary function) as follows. \\
\\
Guess that $y_2(x)=v(x)y_1(x)$ for some $v(x)$. This gives
\begin{align*}
y_2&=vy_1 \\
y_2'&=vy_1'+v'y_1 \\
y_2'&=vy_1''+2v'y_1'+v''y_1.
\end{align*}
We require $y_2$ to be a solution of the equation so (remembering that $y_1$ is also a solution):
\begin{align*}
0&= (vy_1''+2v'y_1'+v''y_1)+p(vy_1' +v'y_1)+q(vy_1) \\
&= v''y_1+v'(2y_1'+py_1)+v(y_1''+py_1'+qy_1) \\
&=v''y_1+v'(2y_1'+py_1)+v(0) \\
&=v''y_1+v'(2y_1'+py_1) \\
\implies \frac{(v')'}{v'}&=-\left( \frac{2y_1'}{y_1}+p \right) \\
\implies \ln v' &=-\left( 2\ln y_1+\int p(x)\mathrm{d} x \right) \\
\implies  v' &=\exp\left[-\left( 2\ln y_1+\int p(x)\mathrm{d} x \right)\right] \\
&=\frac{\exp\left(-\int p(x)\mathrm{d} x \right)}{y_1(x)^2}\\
\implies v&= \int \left[ \frac{\exp\left(-\int p(x)\mathrm{d} x \right)}{y_1(x)^2}\right] \mathrm{d} x \\
\implies y_2(x)&= y_1(x)\int \left[ \frac{\exp\left(-\int p(x)\mathrm{d} x \right)}{y_1(x)^2}\right] \mathrm{d} x \\
\end{align*}
Observe that the small integral gives an arbitrary factor and the large integral adds an arbitrary multiple of $y_1(x)$ so if we keep both degress of freedom, then we have actually found the complete solution.\\
\\
e.g. in the case of a degenerate solution of an equation with constant coefficients, $y_1=e^{kx}$, $p(x)=-2k$ so
\begin{align*}
y&=e^{kx}\int \left[ \frac{\exp\left(-\int (-2k)\mathrm{d} x \right)}{(e^{kx})^2}\right] \mathrm{d} x \\
&=e^{kx}\int \left[ \frac{\alpha e^{2kx}}{e^{2kx}}\right] \mathrm{d} x \\
&=Ae^{kx}\left(x+\frac{\beta}{\alpha}\right) \\
&=e^{kx}\left(\alpha x+\beta\right)
\end{align*}

Q: Phase space
A: Given an $n$th order differential equation, the solution vector \[
\vec{Y}(x)=\left(
\begin{array}{c}
y(x)\\
y'(x)\\
\vdots\\
y^{(n-1)}(x)
\end{array}
\right)
\]
contains all the information necessary to predict the evolution of the system. As $x$ varies, $\vec{Y}(x)$ plots a trajectory through phase space.

Q: The Wronskian
A: The \emph{fundamental matrix} is the matrix $\vec{Y_1},\vec{Y_2}\ldots\vec{Y_n}$ where $\vec{Y_i}$ is the $i$th solution vector. The determinant of this matrix is called the Wronskian.

Q: Abel's theorem
A: \textbf{Abel's Theorem.} Given a 2nd order HLODE \[
y''+p(x)y'+q(x)y=0
\]
where $p$ or $q$ are continuous, then either $W\equiv 0$ or $\forall x: W\neq 0$ .
\begin{proof}
Let $y_1$ and $y_2$ be solutions of the equation such that
\begin{align}
y_1''+p(x)y_1'+q(x)y_1&=0 \label{1}\\
y_2''+p(x)y_2'+q(x)y_2&=0 \label{2}
\end{align}
Thus $y_1\times$\eqref{2}$-y_2\times$\eqref{1} gives
\begin{align*}
0&=y_1 (y_2''+p(x)y_2'+q(x)y_2) -y_2(y_1''+p(x)y_1'+q(x)y_1) \\
&=y_1 y_2''+p(x)y_1y_2' -y_2y_1''-p(x)y_2y_1' \\
&=y_1 y_2''-y_2y_1''+p(y_1y_2' -y_2y_1' )\\
&=W'+pW
\end{align*}
by recognition of $W'=(y_1y_2'-y_2y_1')'=y_1y_2''+y_1'y_2'-y_2'y_1'-y_2y_1''=y_1 y_2''-y_2y_1''$. Thus we have Abel's identity: \[
W=W_0\exp\left[ -\int_{x_0}^xp(u)\mathrm{d}u\right]
\]
where $W_0=W(x_0)=y_1(x_0)y_2'(x_0)-y_2(x_0)y_1'(x_0)$. Thus either $W_0=0\equiv W$ or $W_0\neq 0$ so $\forall x: W\neq 0$.
\end{proof}
In fact, Abel's Theorem holds for $n$th order differentiation equations as they can be expressed as $\vec{Y}'+A(x)\vec{Y}=0$ such that $W'+Tr[A(x)]W=0$ and $W=W_0\exp\left[ -\int_{x_0}^xTr[A(u)]\mathrm{d}u\right]$.

Q: When is the Wronskian constant?
A: Abel's identity states \[
W=W_0\exp\left[ -\int_{x_0}^xp(u)\mathrm{d}u\right]
\]
so $W$ is constant if $W_0=0\equiv W$ or $-\displaystyle\int_{x_0}^xp(u)=c \implies p(u)=0$.

Q: Using the Wronskian to find a second independent solution
A: \begin{align*}
y_1y_2'-y_2y_1'&=W_0\exp\left[-\int_0^xp(u)\mathrm{d}u\right] \\
\implies y_2'-\frac{y_1'(x)}{y_1(x)}y_2&=\frac{W_0\exp\left[-\int_0^xp(u)\mathrm{d}u\right]}{y_1(x)} \\
\implies \frac{y_2'}{y_1}-\frac{y_1'(x)}{y_1^2(x)}y_2&=\frac{W_0\exp\left[-\int_0^xp(u)\mathrm{d}u\right]}{y_1^2(x)} \\
&=\frac{\mathrm{d}}{\mathrm{d}x}\left(\frac{y_2}{y_1}\right) \\
\implies \frac{y_2}{y_1} &= \int \left( \frac{W_0\exp\left[-\int_0^xp(u)\mathrm{d}u\right]}{y_1^2(x)} \right) \mathrm{d}x \\
\implies y &= y_1\int \left( \frac{W_0\exp\left[-\int_0^xp(u)\mathrm{d}u\right]}{y_1^2(x)} \right) \mathrm{d}x \\
\implies y &= y_1\int \left( \frac{\exp\left[-\int_0^xp(u)\mathrm{d}u\right]}{y_1^2(x)} \right) \mathrm{d}x \\
\end{align*}
This is a faster route to the same result achieved by reduction of order.

Q: Magical detuning technique for SHM with resonance
A: Given $\ddot y + \omega_0^2 y = \sin \omega_0 t$, we instead consider the detuned equation  \[
\ddot y + \omega_0^2 y = \sin \omega t \quad (\omega \not= \omega_0 )
\]
Try
\[
  y_p = C(\sin \omega t - \sin \omega_0 t).
\]
because we know the famous $\omega_0^2-\omega^2$ denominator is coming and we want the numerator to be going to zero as $\omega\to\omega_0$ as well, so we subtract a multiple of the complementary function. \\
\\
We have
\[
  \ddot y_p = C(-\omega^2 \sin \omega t + \omega_0^2 \sin\omega_0 t).
\]
Substituting into the differential equation, we have $C(\omega_0^2 - \omega^2) = 1$. Then
\[
  y_p = \frac{\sin \omega t - \sin \omega_0t}{\omega_0^2 - \omega^2}.
\]
Using the factor formula, this becomes
\[
  y_p = \frac{2}{\omega_0^2 - \omega^2} \cos \left(\frac{\omega_0 + \omega}{2}t\right) \sin \left(\frac{\omega - \omega_0}{2} t\right)
\]
We let $\omega_0 - \omega = \Delta \omega$. Then
\begin{align*}
  y_p &= \frac{-2}{(2\omega + \Delta \omega)\Delta \omega}\cos \left[\left(\omega + \frac{\Delta \omega}{2}\right)t\right] \sin \left(\frac{\Delta \omega}{2}t\right) \\
&= -\frac{\cos \left[\left(\omega + \frac{\Delta \omega}{2}\right)t\right]}{(2\omega + \Delta \omega)} \frac{\sin \left(\frac{\Delta \omega}{2}t\right)}{\frac{\Delta \omega}{2}} \\
\implies  \lim_{\omega \to\omega_0,\Delta\omega\to 0} y_p &=-\frac{\cos (\omega_0 t)}{2\omega_0 } \frac{ \frac{\Delta \omega}{2}t}{\frac{\Delta \omega}{2}} \\
&=-\frac{t}{2\omega_0 } \cos (\omega_0 t)\\
\end{align*}
We essentially have an $\omega_0$ frequency sinusoidal wave `beating ' in amplitude with the tiny frequency $\frac{\Delta \omega}{2}$. Thus, as this vanishes, the period of the beats diverges such that the amplitude just keeps increasing.

Q: Variation of parameters
A: Let $y_1(x)$ and $y_2(x)$ be LI CFs of the ODE $y'' + p(x)y' + q(x)y = f(x)$. We take (arbitrarily) the solution vectors $\mathbf{Y}_1 = (y_1, y_1')$ and $\mathbf{Y}_2 = (y_2, y_2')$ to form a basis of the solution space such that
\[
  \mathbf{Y}_p(x) = u(x)\mathbf{Y}_1(x) + v(x)\mathbf{Y}_2(x).
\]
Considering each row explictly:
\begin{align*}
  y_p &= uy_1 + vy_2 \\
  y_p' &= uy_1' + vy_2' \\
  \implies y''_p &= (uy_1'' + u'y_1') + (vy_2'' + v'y_2')
\end{align*}
Substituting into the ODE gives
\begin{align*}
f&=(uy_1'' + u'y_1') + (vy_2'' + v'y_2')  + p( uy_1' + vy_2') + q( uy_1 + vy_2) \\
&=u (y_1'' +py_1+qy_1) + v (y_2''+py_2+qy_2) + y_1' u' + y_2'v'\\
&=u (0) + v (0) + y_1' u' + y_2'v'\\
\implies f&= y_1' u' + y_2'v' \tag{1}\\
\end{align*}
Now we differentiate the first row and compare to the second row:
\begin{align*}
 y_p' &=u'y_1+uy_1'+v'y_2+vy_2' \\
&= uy_1' + vy_2' \\
\implies 0&= u'y_1+v'y_2  \tag {2}
\end{align*}
Thus we have two simultaneous equations for $u'$ and $v'$, which we should be able to solve by matrix methods:
\[
  \begin{pmatrix}
    y_1 & y_2\\
    y_1' & y_2'
  \end{pmatrix}
  \begin{pmatrix}
    u'\\
    v'
  \end{pmatrix}
  =
  \begin{pmatrix}
    0\\
    f
  \end{pmatrix}
\]
Inverting the left matrix (which we can do for LI $y_1,y_2$ by Abel's Theorem), we have
\[
  \begin{pmatrix}
    u'\\
    v'
  \end{pmatrix} = \frac{1}{W}
  \begin{pmatrix}
    y_2' & -y_2\\
    -y_1' & y_1
  \end{pmatrix}
  \begin{pmatrix}
    0\\f
  \end{pmatrix}
\]
So $u' = -\frac{y_2}{W}f$ and $v' = \frac{y_1}{W}f$. Therefore\[
y_p(x)=y_2(x)\int^x \frac{y_1(\xi)f(\xi)}{W(\xi)} \mathrm{d}\xi-y_1(x)\int^x \frac{y_2(\xi)f(\xi)}{W(\xi)} \mathrm{d}\xi
\]
Note finally that the constants are arbitrary as they add multiples of the CFs to the PI.

Q: Solution of equidimensional equation
A: We are given \[
    ax^2y'' + bxy' + cy = 0
  \]
where $a, b, c$ are constants.
\subsubsection*{Solving by eigenfunctions}
Note that $y = x^k$ is an eigenfunction of $x\frac{\mathrm{d}}{\mathrm{d} x}$ which gives
\begin{align*}
y' &= kx^{k - 1} \\
\implies xy' &= kx^k \\
\implies y'' &= k(k - 1)x^{k - 2} \\
\implies x^2y'' &= k(k - 1)x^k
\end{align*}
Substituting in, we have
\[
  ak(k - 1) + bk + c = 0,
\]
which we can solve, in general, to give two roots $k_1$ and $k_2$, and $y_c = Ax^{k_1} + Bx^{k_2}$.
\subsubsection*{Solving by substitution}
Alternatively, we can make a substitution $z = \ln x$. This gives
\begin{align*}
\frac{\mathrm{d}y}{\mathrm{d} x}&=\frac{\mathrm{d}y}{\mathrm{d} z}\frac{\mathrm{d}z}{\mathrm{d} x}\\
&=\frac{1}{x}\frac{\mathrm{d}y}{\mathrm{d} z}\\
\implies x\frac{\mathrm{d}y}{\mathrm{d} x}&=\frac{\mathrm{d}y}{\mathrm{d} z} \\
\implies \frac{\mathrm{d}y}{\mathrm{d} x}+x\frac{\mathrm{d}^2y}{\mathrm{d} x^2}&=\frac{\mathrm{d}^2y}{\mathrm{d} z^2}\frac{\mathrm{d}z}{\mathrm{d} x} \\
\implies \frac{1}{x}\frac{\mathrm{d}y}{\mathrm{d} z} +x\frac{\mathrm{d}^2y}{\mathrm{d} x^2} &=\frac{1}{x}\frac{\mathrm{d}^2y}{\mathrm{d} z^2} \\
\implies  x^2\frac{\mathrm{d}^2y}{\mathrm{d} x^2} &=\frac{\mathrm{d}^2y}{\mathrm{d} z^2} - \frac{\mathrm{d}y}{\mathrm{d} z}
\end{align*}
Substituting in, we have
\begin{align*}
  a \left(\frac{\mathrm{d}^2 y}{\mathrm{d} z^2}-\frac{\mathrm{d}y}{\mathrm{d} z}\right)  + b\frac{\mathrm{d}y}{\mathrm{d} z} + cy &= 0 \\
\implies a \frac{\mathrm{d}^2 y}{\mathrm{d} z^2} + (b-a) \frac{\mathrm{d}y}{\mathrm{d} z} + cy &= 0
\end{align*}
This has solution $y = Ae^{k_1z} + Be^{k_2z} = Ax^{k_1} + Bx^{k_2}$. This method makes it easy to see what happens when there is a degenerate eigenvalue $k=k_1=k_2$:
\begin{align*}
y &= (A+Bz) e^{kz} \\
&= (A+B\ln x) x^k \\
\end{align*}
Similarly, if there is resonance, then there is a particular integral of the form $x^{k_1}\ln x$ or $x^{k_2}\ln x$.

Q: Solving difference equations
A: The eigenfunction of the differential operator $\mathcal{D}[y_n]=y_{n+1}$ is $y_n=\lambda^n$. Thus given \[
ay_{n+2}+by_{n+1}+cy_n=f_n
\]
we solve the homogeneous equation with the ansatz $y_n=\lambda^n$:
\begin{align*}
a\lambda^{n+2}+b\lambda^{n+1}+c\lambda^n&=0 \\
\implies a\lambda^2+b\lambda+c&=0
\end{align*}
which is a characteristic quadratc equation with two roots. Therefore, the general complementary function is \[
y_n^{(c)}=\begin{cases}
  A\lambda_1^n+b\lambda_2^n  &\textrm{if } \lambda_1\neq\lambda_2 \\
 (A+Bn)\lambda^n &\textrm{if }\lambda=\lambda_1=\lambda_2
\end{cases}
\]
Finally, the particular `integral' is \[
y_n^{(p)}=\begin{cases}
A\lambda^n &\textrm{if } f_n=\lambda^n,\lambda\neq\lambda_1,\lambda_2 \\
An\lambda_1^n &\textrm{if }f_n=\lambda_1^n \\
Bn\lambda_2^n &\textrm{if }f_n=\lambda_2^n \\
\sum_{r=0}^p a_r n^r &\textrm{if }f_n=n^p
\end{cases}
\]

Q: Forced SHM with damping
A: \[
M\ddot{x}=F-kx-l\dot{x}
\]
As $kx$ is a force, we know it has dimesions $\frac{\textrm{mass}\times\textrm{distance}}{ \textrm{time}^2}$ so $k$ has dimensions $\frac{\textrm{mass}}{\textrm{time}^2}$. Therefore, a natural time scale to consider is $t_c=\sqrt{\frac{M}{k}}$. Thus let us define a nondimensional time by $t=\tau t_c$ such that
\begin{align*}
\frac{\mathrm{d}}{\mathrm{d}t}&=\frac{1}{t_c}\frac{\mathrm{d}}{\mathrm{d}\tau} \\
&=\sqrt{\frac{k}{M}}\frac{\mathrm{d}}{\mathrm{d}\tau}
\end{align*}
This gives the rescaled equation
\begin{align*}
\frac{k}{M} \left(M \frac{\mathrm{d}^2x}{\mathrm{d}\tau^2} \right)&=F-kx-l \sqrt{\frac{k}{M}}\frac{\mathrm{d}x}{\mathrm{d}\tau}  \\
\implies \frac{\mathrm{d}^2x}{\mathrm{d}\tau^2} +\frac{l}{\sqrt{kM}}\frac{\mathrm{d}x}{\mathrm{d}\tau}+x &=\frac{F}{k} \\
\implies \frac{\mathrm{d}^2x}{\mathrm{d}\tau^2} +2\kappa\frac{\mathrm{d}x}{\mathrm{d}\tau}+x &=f \\
\end{align*}
where $f=\frac{F}{k}$ and $\kappa=\frac{l}{2\sqrt{kM}}$. This gives the characteristic equation
\begin{align*}
\lambda^2+2\kappa\lambda+1=0 \\
\implies (\lambda+\kappa)^2+1-\kappa^2=0 \\
\implies \lambda=-\kappa\pm\sqrt{\kappa^2-1}
\end{align*}
Therefore \[
x=
\begin{cases}
e^{-\kappa\tau} \left( A\sin\sqrt{1-\kappa^2}t+B\cos\sqrt{1-\kappa^2}t \right) &\textrm{ if $\kappa<1$ (underdamping)} \\
(A+B\tau) e^{-\kappa\tau} &\textrm{ if $\kappa=1$ (critical damping)} \\
Ae^{\lambda_1t}+Be^{\lambda_2t};\lambda_{1,2}=-\kappa \pm\sqrt{\kappa^2-1}&\textrm{ if $\kappa>1$ (overdamping)}
\end{cases}
\]
In underdamping, the system oscillates with period $2\pi/\sqrt{1-\kappa^2}$ and the amplitude decays exponentially in time $O(1/\kappa)$. In overdamping, note that $\lambda_{1,2}<0$ but $\lambda_2$ is more negative, so the second term dominates initially. Therefore, if $B<0$, as the second negative term exponentially decays, there is a sudden increase in the response for a time $O(1/\lambda_2)$. Then the first term dominates and decays in time $O(1/\lambda_1)$. \\
\\
Critical damping is the limit of underdamping or overdamping. The period of the underdamping goes to infinity leaving just the first wave cycle never quite finishing. The graph looks the same as overdamping if you remember that $\lambda_{1,2}\to\kappa=1$  because initial the linear factor dominates and then there is an exponential decay. \\
\\
If there is forcing, then the complementary functions dominate initially (transient response) but then decay, leaving the particular integral to determine the ultimate behaviour (asymptotic response.

Q: Dynamics of a bouncing ball
A: Resolving forces \[
m\ddot{x}=F(t)-mg
\]
where $F(t)$ is an impulsive force, i.e. it acts instantaneously to push the ball back up. Thus, integrating:
\begin{align*}
\int_{T-\epsilon}^{T+\epsilon} m \frac{\mathrm{d}^2x}{\mathrm{d}t^2} \mathrm{d}t= \int_{T-\epsilon}^{T+\epsilon} F(t) \mathrm{d}t - \int_{T-\epsilon}^{T+\epsilon} mg \mathrm{d}t \\
\left[ m\frac{\mathrm{d}x}{\mathrm{d}t}\right]^{T+\epsilon}_{T-\epsilon}=I-2mg\epsilon \\
\left[ m\frac{\mathrm{d}x}{\mathrm{d}t}\right]^{T+}_{T-}=I \\
\end{align*}
Thus the only property of the force which must be known to determine the evolution of the system (in particular the change in momentum) is the impulse \[
I= \int_{T-\epsilon}^{T+\epsilon} F(t) \mathrm{d}t
\]
This can be conveniently represented by the equation \[
m\ddot{x}=I\delta(t-T)-mg
\]
where the impulse acts at time $T$.

Q: Dirac $\delta$-function
A: Consider a family of functions $\mathcal{D}(t,\epsilon)$ such that
\begin{enumerate}
\item \[ \forall t\neq 0: \lim_{\epsilon\to 0} \mathcal{D}(t,\epsilon) =0 \]
\item \[ \lim_{\epsilon\to 0} \int_{-\infty}^\infty \mathcal{D}(t,\epsilon) =1 \]
\end{enumerate}
Now let \[
\delta(x)=\lim_{\epsilon\to 0} \mathcal{D}(t,\epsilon)
\]
which is not a defined function, but instead is a shorthand for its integral properties.
\begin{enumerate}
\item $ \forall x\neq 0: \delta(x)=0$
\item For any test function $g(x)$ continuous at $x=x_0$: \[
\int_a^b g(x)\delta(x-x_0) =
\begin{cases}
g(x_0) \quad &\textrm{if } a<x_0<b \\
0 \quad &\textrm{if } x_0<a \textrm{ or }x_0<b
\end{cases}
\]
\item $$ \int_{-\infty}^\infty \delta(x)=1$$
\end{enumerate}

Q: $\delta$-function forcing
A: Given
\[
ay''+by'+cy=I\delta(x-x_0)
\]
we first solve
\[
ay''+by'+cy'=0
\]
for $x<x_0$ and then the same equation for $x>x_0$. Now, integrating from $x_0^-$ to $x_0^+$:
\begin{align*}
a[y']_{x_0^-}^{x_0^+}+b[y]_{x_0^-}^{x_0^+}+c\int_{x_0^-}^{x_0^+} y \mathrm{d}x=I \\
\implies a[y']_{x_0^-}^{x_0^+}=I
\end{align*}
as $y$ is continuous (else $y''$ would be too discontinuous for the equation to hold). Thus we use the jump conditions
\begin{align*}
a[y']_{x_0^-}^{x_0^+}=I \\
[y]_{x_0^-}^{x_0^+}=0
\end{align*}
to solve for the last two initial conditions.

Q: The Heaviside Step function
A: We define \[
H(x)=\int_{-\infty}^x\delta(t)\mathrm{d}t
\]
which gives  \[
H(x)= \begin{cases}
0 \quad &\textrm{for }  x<0 \\
1 \quad &\textrm{for } x>0  \\
\textrm{undefined} \quad &\textrm{at } x=0
\end{cases}
\]
and \[
\frac{\mathrm{d}H}{\mathrm{d}x}=\delta(x)
\]

Q: Heaviside step function forcing
A: Given \[
by'+cy=IH(x-x_0),
\]
if we consider the change from $x_0^-$ to $x_0^+$:
\begin{align*}
b[y']_{x_0^+}^{x_0^-}+c[y]_{x_0^+}^{x_0^-}=I \\
\implies b[y']_{x_0^+}^{x_0^-}=I
\end{align*}
This gives the jump conditions
\begin{align*}
b[y']_{x_0^+}^{x_0^-}=I \\
[y]_{x_0^+}^{x_0^-}=0
\end{align*}

Q: Classification of singular points
A: We will consider equations of the form
\[
  p(x) y'' + q(x) y' + r(x) y = 0.
\]
   If $\frac{q}{p}$ and $\frac{r}{p}$ have Taylor series about $x_0$, then the point $x = x_0$ is an \emph{ordinary point}. Otherwise, $x_0$ is a \emph{singular point}. \\
\\
  If $x_0$ is a singular point but the equation can be written as
  \[
    P(x)(x - x_0)^2y'' + Q(x)(x - x_0)y' + R(x)y = 0,
  \]
  where $\frac{Q}{P}$ and $\frac{R}{P}$ have Taylor series about $x_0$, then $x_0$ is a \emph{regular singular point} (aka removable singularity). That is, if multiplying $\frac{q}{p}$ by $x-x_0$ and $\frac{r}{p}$ by $(x-x_0)^2 $ means that they now have Taylor series expansions, then the singular point is ordinary. Otherwise, $x_0$ is an irregular singular point (aka essential singularity).

Q: Series solution of a second order LHODE about an ordinary point
A: Suppose we have \[
  p(x) y'' + q(x) y' + r(x) y = 0.
\]
for which $x=x_0$ is an ordinary point. Try
\begin{align*}
y&=\sum_{n=0}^\infty a_n (x-x_0)^n \\
\implies (x-x_0)y'&=\sum_{n=1}^\infty na_n(x-x_0)^n \\
\implies (x-x_0)^2y''&=\sum_{n=2}^\infty n(n-1)a_n(x-x_0)^n \\
\end{align*}
Thus
 \begin{align*}
0&=p(x-x_0)^2y''+q(x-x_0) [(x-x_0)y']+r(x-x_0)^2y \\
&=p\sum_{n=2}^\infty n(n-1)a_n(x-x_0)^n + q(x-x_0)\sum_{n=1}^\infty na_n(x-x_0)^n +r(x-x_0)^2 \sum_{n=0}^\infty  a_n(x-x_0)^n\\
&=\sum_{n=2}^\infty n(n-1)a_np(x-x_0)^n + \sum_{n=1}^\infty na_n[q(x-x_0)] (x-x_0)^n + \sum_{n=0}^\infty  a_n[r(x-x_0)^2 ](x-x_0)^n \\
&=\sum_{n=2}^\infty \left[ n(n-1)a_np +(n-1)a_{n-1} q + a_{n-2}r \right] (x-x_0)^n
%\sum_{n=2}^\infty
\end{align*}
Equating coefficients of $(x-x_0)^n$ gives a recurrence relation between the $a_n$ (at best for $n\geqslant 2$, but if $p(x),q(x),r(x)$ have some higher order terms then you may have to consider some small values of $n$ separately). The coefficients $a_0$ and $a_1$ are arbitrary and the other coefficients must be determined in terms of those two using the recurrence relation many (usually $n$) times. Finally, you may be able to recognise the result as a Taylor series.

Q: Series solution of a second order LHODE about a regular singular point
A: Suppose we have \[
  p(x) y'' + q(x) y' + r(x) y = 0.
\]
for which $x=x_0$ is an regular singular point. By Fuchs' Theorem we are guaranteed at least one solution of the form \[
y=\sum_{n=0}^\infty a_n (x-x_0)^{n+\sigma}
\]
where $a_0\neq 0$ by construction and $\sigma\in\mathbb{R}$.
Try
\begin{align*}
y&=\sum_{n=0}^\infty a_n (x-x_0)^{n+\sigma} \\
\implies (x-x_0)y'&=\sum_{n=1}^\infty (n+\sigma)a_n(x-x_0)^{n+\sigma} \\
\implies (x-x_0)^2y''&=\sum_{n=2}^\infty (n+\sigma)(n+\sigma-1)a_n(x-x_0)^{n+\sigma} \\
\end{align*}
Thus
 \begin{align*}
0&=p(x-x_0)^2y''+q(x-x_0) [(x-x_0)y']+r(x-x_0)^2y \\
&=p\sum_{n=2}^\infty (n+\sigma)(n+\sigma-1)a_{n+\sigma}(x-x_0)^{n+\sigma} + q(x-x_0)\sum_{n=1}^\infty na_n(x-x_0)^n +r(x-x_0)^2 \sum_{n=0}^\infty  a_n(x-x_0)^n
\end{align*}
Now we choose $n$ such that we get the coefficient of $a_0$ and consider the smallest power of $x$. Recalling that we construct $a_0$ to be non-zero, this gives a quadratic equation for $\sigma$ known as the indicial equation. \\
\\
Case 1: $\sigma_1-\sigma_2$ is a non-zero non-integer. We are guaranteed that these two values of $\sigma$ lead to 2 LI solutions of the ODE. For each value of $\sigma$, equating coefficients of $(x-x_0)^{n+\sigma}$ gives a recurrence relation between the $a_n$ (at best for $n\geqslant 2$, depending on $p(x),q(x),r(x)$). The coefficient $a_0$ is arbitrary. The coefficient $a_1$ is found by letting $n=1$ and considering the smallest power of $x$. The other coefficients must be determined in terms of $a_0$ using the recurrence relation many (usually $n$) times. \\
\\
Case 2: $\sigma_1-\sigma_2$ is a non-zero integer. Let $\operatorname{Re}(\sigma_1)>\operatorname{Re}(\sigma_2)$. Then $y_1 =\sum_{n=0}^\infty a_n (x-x_0)^{n+\sigma_1}$ and $y_2=\sum_{n=0}^\infty b_n (x-x_0)^{n+\sigma_2}+c\ln(x-x_0)\sum_{n=0}^\infty a_n (x-x_0)^{n+\sigma_1}$. (Constructing the logarithmic solution requires reduction of order.) $c$ has to be non-zero only if the smaller root has its first coefficient zero such that it simply recovers the first series solution. If there is any LI part of the second series solution then that is sufficient. \\
\\
Case 3: $\sigma_1=\sigma_2$. Same as case $2$ but we are guaranteed $c\neq 0$.

Q: Proof that second order LHODEs have at most two LI solutions
A: Suppose $y_1,y_2,y_3$ are solutions of \[
y''+py'+qy=0
\]
Therefore:
\begin{align*}
\left|\begin{array}{ccc}
y_1 & y_2 & y_3 \\
y_1' & y_2' & y_3' \\
y_1'' & y_2'' & y_3''
\end{array} \right| \\
\end{align*}
is
\begin{align*}
&\left|\begin{array}{ccc}
y_1 & y_2 & y_3 \\
y_1' & y_2' & y_3' \\
-(py_1'+qy_1)' & -(py_2'+qy_2) & -(py_3'+qy_3)
\end{array} \right| \\
&=y_1(-y_2'(py_3'+qy_3)+y_3'(py_2'+qy_2) -y_2 (-y_1'(py_3'+qy_3)+y_3'(py_1'+qy_1)) +y_3 (-y_1'(py_2'+qy_2)+y_2'(py_1'+qy_1) )\\
&=0
\end{align*}
Thus the solutions are linearly dependent. \qed

Q: Directional derivative
A: From the chain rule
\begin{align*}
\mathrm{d}f&=\frac{\partial f}{\partial x} \mathrm{d}x+\frac{\partial f}{\partial y} \mathrm{d}y  \\
&= (\mathrm{d}x,\mathrm{d}y) \cdot \left( \frac{\partial f}{\partial x}, \frac{\partial f}{\partial y} \right) \\
&=\vec{ds}\cdot\nabla f \\
&=(\mathrm{d}s \ \hat{s})\cdot\nabla f \\
&=\mathrm{d}s(\hat{s}\cdot\nabla f ) \\
&=\mathrm{d}s|\nabla f|\cos\theta
\end{align*}
This prompts the definition of the directional derivative of $f$ with respect to $s$ \[
\frac{\mathrm{d}f}{\mathrm{d}s}=\hat{s}\cdot\nabla f
\]

Q: Gradient: definition and properties
A: We define\[
\nabla f =(\mathrm{d}x,\mathrm{d}y).
\]
As \[
\frac{\mathrm{d}f}{\mathrm{d}s}=|\nabla f|\cos\theta
\]
we can see that as $\theta$ varies, the directional derivative is maximised at $\theta=0$ such that \[
\max_\theta \left| \frac{\mathrm{d}f}{\mathrm{d}s}\right|=|\nabla f|.
\]
Therefore, $\nabla f$ is a vector with the following properties:
\begin{itemize}
\item \emph{magnitude} equal to the maximum rate of change of $f$ with distance
\item pointing in the \emph{direction} of the greatest increase of $f$
\item If $\hat{s}$ is a contour then \[ \frac{\mathrm{d}f}{\mathrm{d}s} =\hat{s}\cdot\nabla f =0\] so $\nabla f$ is everywhere orthogonal to \emph{contours}
\end{itemize}

Q: Condition for stationary points of a multivariate function
A: \begin{align*}
\forall\hat{s}: \ &\hat{s}\cdot\nabla f=0 \\
\implies &\nabla f=\vec{0} \\
\implies &\frac{\partial f}{\partial x}=\frac{\partial f}{\partial y}=0
\end{align*}

Q: Multivariate Taylor series
A: Suppose we have a function $f(x, y)$ and a point $\mathbf{x}_0$. Now consider a finite displacement $\delta s$ along a straight line in the $x,y$ plane. From the directional derivative formula
\[
  \delta s\frac{\mathrm{d} f}{\mathrm{d} s} = \delta \mathbf{s}\cdot \nabla f
\]
The Taylor series along the line is
\begin{align*}
  f(\vec{s}) &= f(\vec{s_0} + \vec{\delta s})\\
  &=f(\vec{s_0})  + \delta s\frac{\mathrm{d} f}{\mathrm{d} s} + \frac{1}{2}(\delta s)^2\frac{\mathrm{d} ^2 f}{\mathrm{d} s^2}\\
  &= f(\vec{s_0}) + \delta \mathbf{s} \cdot \nabla f + \frac{1}{2}(\delta s)^2 (\mathbf{\hat{s}}\cdot \nabla)(\mathbf{\hat{s}}\cdot \nabla)f.
\end{align*}
Where
\begin{align*}
  \delta \mathbf{s}\cdot \nabla f &= (\delta x)\frac{\partial f}{\partial x} + (\delta y)\frac{\partial f}{\partial y}\\
  &= (x - x_0)\frac{\partial f}{\partial x} + (y - y_0)\frac{\partial f}{\partial y}
\end{align*}
and
\begin{align*}
  \delta s^2 (\mathbf{\hat{s}}\cdot \nabla)(\mathbf{\hat{s}}\cdot \nabla)f &= (\delta \mathbf{s}\cdot \nabla)(\delta \mathbf{s}\cdot \nabla) f\\
  &= \left[\delta x\frac{\partial }{\partial x} + \delta y\frac{\partial}{\partial y}\right]\left[\delta x\frac{\partial f}{\partial x} + \delta y\frac{\partial f}{\partial y}\right]\\
  &= (\delta x)^2 \frac{\partial^2 f}{\partial x^2} + 2\delta x\delta y\frac{\partial^2 f}{\partial x\partial y} + (\delta y)^2 \frac{\partial^2 f}{\partial y^2}\\
  &=
  \begin{pmatrix}
    \delta x& \delta y
  \end{pmatrix}
  \begin{pmatrix}
    f_{xx}&f_{xy}\\
    f_{yx}&f_{yy}
  \end{pmatrix}
  \begin{pmatrix}
    \delta x\\\delta y
  \end{pmatrix}
\end{align*}
  The \emph{Hessian matrix} is the matrix
  \[
    \mathsf{H}=\nabla \nabla f =
    \begin{pmatrix}
      f_{xx}&f_{xy}\\
      f_{yx}&f_{yy}
    \end{pmatrix}
  \]
In conclusion, we have
\begin{align*}
  f(x, y) &= f(x_0, y_0) + (x - x_0)f_x + (y - y_0)f_y \\
  &+ \frac{1}{2}[(x - x_0)^2 f_{xx} + 2(x - x_0)(y - y_0)f_{xy} + (y - y_0)^2 f_{yy}]
\end{align*}
In general, the coordinate-free form is
\[
  f(\mathbf{s}) = f(\mathbf{s_0}) + \delta \mathbf{s}\cdot\nabla f(\mathbf{s}_0) + \frac{1}{2}\delta \mathbf{s}\cdot \nabla \nabla f\cdot \delta \mathbf{s}
\]
where the dot in the second term represents a matrix product. Alternatively, in terms of the gradient operator (and real dot products), we can write
\[
  f(\mathbf{s}) = f(\mathbf{s_0}) + \delta \mathbf{s}\cdot \nabla f(\mathbf{s}_0) + \frac{1}{2}[\nabla (\nabla f\cdot \delta \mathbf{s})]\cdot \delta\mathbf{s}
\]

Q: Classifying stationary points of a multivariate function
A: At a stationary point $\mathbf{s}_0$, we know that $\nabla f(\mathbf{s}_0) = 0$. So at a point near the stationary point,
\[
  f(\mathbf{s})\approx f(\mathbf{s}_0) + \frac{1}{2}\delta\mathbf{s}\cdot H\cdot \delta\mathbf{s},
\]
where $H = \nabla\nabla f(\mathbf{s}_0)$ is the Hessian matrix. \\
\\
At a minimum, every point near $\mathbf{s}_0$ has $f(\mathbf{s} )> f(\mathbf{s}_0)$, ie. $\delta \mathbf{s}\cdot H\cdot \delta \mathbf{s} > 0$ for all $\delta \mathbf{s} (\neq 0)$. We say $\delta \mathbf{s}\cdot H\cdot\delta \mathbf{s}$ is \emph{positive definite}.
Similarly, at a maximum, $\delta \mathbf{s}\cdot H\cdot \delta\mathbf{s} < 0$ for all $\delta \mathbf{s} (\neq 0)$. We say $\delta \mathbf{s}\cdot H\cdot \delta \mathbf{s}$ is \emph{negative definite}.
At a saddle, $\delta \mathbf{s} \cdot H\cdot \delta\mathbf{s}$ is indefinite, ie. it can be positive, negative or zero depending on the direction.
This, so far, is not helpful, since we do not have an easy way to know what sign $\delta \mathbf{s}\cdot H\cdot \delta \mathbf{s}$ could be. Now note that $H = \nabla\nabla f$ is symmetric (because $f_{xy} = f_{yx}$). So $H$ can be diagonalized. With respect to the axes in which $H$ is diagonal (principal axes), we have
\begin{align*}
  \delta\mathbf{s}\cdot H\cdot \delta\mathbf{s} &= (\delta x_1, \delta x_2, \cdots, \delta x_n)
  \begin{pmatrix}
    \lambda_1\\
    &\lambda_2\\
    &&\ddots\\
    &&&\lambda_n
  \end{pmatrix}
  \begin{pmatrix}
    \delta x\\\delta y\\\vdots\\\delta z
  \end{pmatrix}\\
  &= \lambda_1(\delta x_1)^2 + \lambda_2 (\delta x_2)^2 + \cdots + \lambda_n(\delta x_n)^2
\end{align*}
where $\lambda_1, \lambda_2, \cdots \lambda_n$ are the eigenvalues of $H$.
So for $\delta \mathbf{s}\cdot H\cdot \delta \mathbf{s}$ to be positive-definite, we need $\lambda_i > 0$ for all $i$. Similarly, it is negative-definite iff $\lambda_i < 0$ for all $i$. If eigenvalues have mixed sign, then it is a saddle point.
Finally, if there is at least one zero eigenvalue and the rest have the same sign, then we need to go to higher orders. \\
Apart from finding eigenvalues, another way to determine the definiteness  is using the \emph{signature}. \\
\\
  The \emph{signature} of $H$ is the pattern of the signs of the subdeterminants:
  \[
    \underbrace{f_{x_1x_1}}_{|H_1|},
    \underbrace{
      \begin{vmatrix}
        f_{x_1x_1} & f_{x_1x_2}\\
        f_{x_2x_1} & f_{x_2x_2}
      \end{vmatrix}}_{|H_2|},\cdots,
    \underbrace{\begin{vmatrix}
      f_{x_1x_1} & f_{x_1x_2} & \cdots & f_{x_1x_n}\\
      f_{x_2x_1} & f_{x_2x_2} & \cdots & f_{x_2x_n}\\
      \vdots & \vdots & \ddots & \vdots\\
      f_{x_nx_1} & f_{x_nx_2} & \cdots & f_{x_nx_n}
    \end{vmatrix}}_{|H_n| = |H|}
  \]
  $H$ is positive definite if and only if the signature is $+, +, \cdots, +$ (as $H$ is a minimum iff. it is a minimum in every subspace). $H$ is negative definite if and only if the signature is $-, +, \cdots, (-1)^n$ (as adding a negative row multiplies the determinant by -1). Otherwise, $H$ is indefinite. \\
\\
If any of the determinants are zero and they have the same sign, then once again higher orders are required.

Q: Contours in the neighbourhood of a stationary point
A: Consider $H$ in 2 dimensions, and axes in which $H$ is diagonal. So $H =
\begin{pmatrix}
  \lambda_1 & 0\\
  0 & \lambda_2
\end{pmatrix}$. Write $\mathbf{s} - \mathbf{s}_0 = (X, Y)$. \\
\\
Then near $\mathbf{s}_0$:
\begin{align*}
f(\vec{s})&=f(\vec{s}_0)+ \mathbf{\delta s}\cdot H \cdot\mathbf{\delta s} \\
&=\textrm{constant}  \\
\implies \mathbf{\delta s}\cdot H \cdot\mathbf{\delta s} &= \textrm{constant}  \\
\implies \lambda_1 X^2 + \lambda_2 Y^2 &= \textrm{constant}
\end{align*}
At a maximum or minimum, $\lambda_1$ and $\lambda_2$ have the same sign. So these contours are locally ellipses. At a saddle point, they have different signs and the contours are locally hyperbolae.

Q: Sketching contour maps
A: \begin{enumerate}
\item Find and classify stationary points
\item Draw ellipses and hyperbolae around stationary points
\item Consider extreme values of $x$ and $y$
\item Consider direction of $\nabla f$ and recall that contours are perpendicular to gradients
\item Picking simple constant values of $f$ and solving algebraically
\end{enumerate}

Q: Equivalence of coupled first order LODEs and second order LODEs
A: Consider two dependent variables $y_1(t), y_2(t)$ related by
\begin{align*}
  \dot y_1 &= ay_1 + by_2 + f_1(t)\\
  \dot y_2 &= cy_1 + dy_2 + f_2(t)
\end{align*}
We can write this in vector notation by
\[
  \begin{pmatrix}
    \dot y_1\\\dot y_2
  \end{pmatrix} =
  \begin{pmatrix}
    a & b\\
    c & d
  \end{pmatrix}
  \begin{pmatrix}
    y_1\\ y_2
  \end{pmatrix}
+
  \begin{pmatrix}
    f_1\\f_2
  \end{pmatrix}
\]
or $\mathbf{\dot Y} = \mathsf{M}\mathbf{Y} + \mathbf{F}$. We can convert this to a higher-order equation by differentiation:
\begin{align*}
  \ddot y_1 &= a\dot y_1 + b\dot y_2 + \dot f_1\\
  &= a\dot y_1 + b(cy_1 + dy_2 + f_2) + \dot f_1\\
  &= a\dot y_1 + bcy_1 + d(\dot y_1 - ay_1 - f_1) + bf_2 + \dot f_1
\end{align*}
so
\[
  \ddot y_1 - (a + d)\dot y_1 + (ad - bc) y_1 = bf_2 - df_1 + \dot f_1
\]
Conversely, if $\ddot y + a\dot y + by = f$, write $y_1 = y$ and $y_2 = \dot y$.
Our system of equations becomes
\begin{align*}
  \dot y_1 &= y_2\\
  \dot y_2 &= f - a y_2 - by_1
\end{align*}
or
\[
 \begin{pmatrix}
  \dot{y_1}\\ \dot{y_2}
\end{pmatrix}=
  \begin{pmatrix}
    0 & 1\\
    -b & -a
  \end{pmatrix}
  \begin{pmatrix}
    y_1\\y_2
  \end{pmatrix}
  \begin{pmatrix}
    0 \\ f
  \end{pmatrix}.
\] \\
In both cases, note that the second order ODE is non-trivially second order iff. the matrix is invertible. \\
\\
In fact any $n$th order differential equation can be written as a system of $N$ first order differential equations.

Q: Solving systems of differential equations
A: Given \[
\dot{\vec{Y}}=\mathsf{M}\vec{Y}+\vec{F}
\]
we first solve the homogeneous equation \[
\dot{\vec{Y}}=\mathsf{M}\vec{Y}
\]
with the complementary function \[
\vec{Y}_c=\vec{v}e^{\lambda t}
\]
where $\vec{v}$ is a constant vector, as this is the eigenfunction of the differential operator on $n$-dimensional vectors. Therefore \[
\lambda\vec{v}=\mathsf{M}\vec{v}
\]
so the general solution is \[
\vec{Y}_c=\sum_{i=1}^nA_i\vec{v_i}e^{\lambda_i t}
\]
where $\lambda_i,\vec{v}_i$ are the $n$ eigenvalue/ eigenvector pairs of the matrix $\mathsf{M}$ and $A_i$ are arbitrary scalars (as any scalar multiple of an eigenvector is also an eigenvector). \\
\\
The phase space trajectory can then be sketched by drawing the eigenvector lines with arrows pointing towards the origin if their eigenvalue is negative (the real part) and arrows pointing away from the origin if their eigenvalue is positive (the real part). These are the lines of actions of the force affecting the trajectory of a point in phase space. In 2 dimensions, this gives the following cases:
\begin{enumerate}
\item $\lambda_1,\lambda_2 \in \mathbb{R}$ and $\lambda_1\lambda_2<0$: a saddle
\item $\lambda_1,\lambda_2 \in \mathbb{R}$ and $\lambda_1>0, \lambda_2>0$: an unstable node
\item $\lambda_1,\lambda_2 \in \mathbb{R}$ and $\lambda_1<0, \lambda_2<0$: a stable node
\item $\lambda_1=\lambda_2*=\lambda$: a spiral (N.B You must use Euler's identity, recall that the eigenvectors come in complex conjugate pairs and split up one of the solutions into its real and imaginary parts. The other solution will thus be the complex conjugate, so (using the fact that the sum and difference of CFs are CFs) we can always use the real and imaginary parts of the first solution as a basis. To deduce whether the spiral is clockwise or anti-clockwise, find the velocity at a point).
\subitem $\operatorname{Re}(\lambda)>0\implies$ unstable spiral
\subitem $\operatorname{Re}(\lambda)<0\implies$ stable spiral
\subitem $\operatorname{Re}(\lambda)=0\implies$ centre (locally elliptical as it surrounds a stationary point)
\end{enumerate}
Of these generic behaviours, all except the centre are said to be structurally stable, in the sense that small changes in the parameters or elements of the key matrix $\mathsf{M}$ will not change the character of the fixed point. The centre is structurally unstable since a small change is highly likely to introduce a real part into the eigenvalues, thus changing the centre into a spiral. \\
\\
You don't have to worry about zero or repeated eigenvalues. These shapes are all centred around a stationary point. As normal, for the particular integral, guess a function of a similar form to $\vec{F}$.

Q: Fixed points of homogeneous coupled differential equations
A: Consider the second-order autonomous system
\begin{align*}
  \dot x &= f(x, y)\\
  \dot y &= g(x, y)
\end{align*}
It can be difficult to solve the equations if $f,g$ are non-linear, but we can learn a lot about phase-space trajectories of these solutions by studying the equilibria and their stability. \\
\\
  An \emph{equilibrium point} is a point in which $\dot x = \dot y = 0$ at $\mathbf{s}_0 = (x_0, y_0)$.
Clearly this occurs when $f(x_0, y_0) = g(x_0, y_0) = 0$. We solve these simultaneously for $x_0, y_0$.
To determine the stability, write $x = x_0 + \xi$, $y = y_0 + \eta$. Then
\begin{align*}
  \dot \xi &= f(x_0 + \xi, y_0 + \eta)\\
  &= f(x_0, y_0)+ \xi \frac{\partial f}{\partial x}(\mathbf{s}_0) + \eta \frac{\partial f}{\partial y}(\mathbf{s}_0) + O(\xi^2, \eta^2)
\end{align*}
So if $\xi, \eta \ll 1$,
\[
  \begin{pmatrix}
    \dot \xi\\\dot \eta
  \end{pmatrix} =
  \begin{pmatrix}
    f_x & f_y\\
    g_x & g_y
  \end{pmatrix}
  \begin{pmatrix}
    \xi\\\eta
  \end{pmatrix}
\]
which is a linear system (so you can use eigenvalues and eigenvectors) for the behaviour near the fixed points.

Q: First order wave equation
A: \[
  \frac{\partial y}{\partial t} = c\frac{\partial y}{\partial x},
\]
with $c$ a constant and $y$ a function of $x$ and $t$. We will later see that solutions correspond to waves travelling in one direction.
We write this as
\[
  \frac{\partial y}{\partial t} - c\frac{\partial y}{\partial x} = 0.
\]
Recall that along a path $x = x(t)$ so that $y = y(x(t), t)$,
\begin{align*}
  \frac{\mathrm{d} y}{\mathrm{d}t} &= \frac{\partial y}{\partial x}\frac{\mathrm{d} x}{\mathrm{d} t} + \frac{\partial y}{\partial t}\\
  &= \frac{\mathrm{d} x}{\mathrm{d} t}\frac{\partial y}{\partial x} + c\frac{\partial y}{\partial x}
\end{align*}
by the chain rule. Now we choose a path along which
\[
  \frac{\mathrm{d} x}{\mathrm{d} t} = -c. \tag{1}
\]
such that the path is a characteristic of the PDE i.e.
\[
  \frac{\mathrm{d} y}{\mathrm{d} t} = 0 \tag{2}
\]
Each path that satisfies (1) can be described by $x = x_0 - ct$, where $x_0$ is a constant. We can write this as $x_0 = x + ct$.
From (2), along each of these paths, $y$ is constant. Therefore, $y=f(x_0)$. Therefore the solution to the wave equation is
\[
  y = f(x + ct),
\]
where $f$ is determined by initial conditions. By differentiating this directly, we can easily check that every function of this form is a solution to the wave equation. \\
\\
If there is forcing $f(t)$, then \[
\frac{\mathrm{d}y}{\mathrm{d}t}=f(t)
\]
so the characteristic is not a contour of $y$. \\
\\
Note that as we move up the time axis, we are simply taking the $t = 0$ solution and translating it to the left. Thus, the wave is moving to the left.

Q: Second-order homogeneous wave equation (``hyperbolic equation'') general solution
A: \[
  \frac{\partial ^2 y}{\partial t^2} = c^2 \frac{\partial^2 y}{\partial x^2}.
\]
Let $y(x, t)$ be the vertical displacement of a horizontal string at the point $x$ at time $t$. The equation says the acceleration at a point is proportional to the curvature of the string at that point.
\begin{align*}
  \frac{\partial ^2 y}{\partial t^2} - c^2 \frac{\partial^2 y}{\partial x^2} &= 0\\
  \left(\frac{\partial}{\partial t} + c\frac{\partial}{\partial x}\right)\left(\frac{\partial}{\partial t} - c\frac{\partial}{\partial x}\right) y &= 0
\end{align*}
Since the operators are commutative, $y = f(x - ct)$ and $y=g(x+ct)$ are both solutions. Since the equation is linear, the general solution is
\[
  y = f(x + ct) + g(x - ct).
\]
This shows that the solution is a superposition of a wave travelling to the left and waves travelling to the right. We can show that this is indeed the most general solution by substituting $\xi = x + ct$ and $\eta = x - ct$.
\begin{align*}
\frac{\partial}{\partial x}&=\frac{\partial\xi}{\partial x}\frac{\partial}{\partial \xi}+\frac{\partial\eta}{\partial x}\frac{\partial}{\partial \eta}  \\
&=\frac{\partial}{\partial \xi}+\frac{\partial}{\partial \eta} \\
\implies c^2 \frac{\partial^2y}{\partial x^2}&=c^2 \left( \frac{\partial}{\partial \xi}+\frac{\partial}{\partial \eta}\right)\left( \frac{\partial}{\partial \xi}+\frac{\partial}{\partial \eta}\right)y \\
&=c^2 \left( \frac{\partial^2}{\partial \xi^2} +\frac{\partial^2}{\partial \xi\partial\eta}+ \frac{\partial^2}{\partial \eta\partial\xi}  +\frac{\partial^2}{\partial \eta^2}\right)y \\
\frac{\partial}{\partial t}&=c\frac{\partial}{\partial \xi}-c\frac{\partial}{\partial \eta}  \\
\implies \frac{\partial^2y}{\partial t^2}&= \left( c\frac{\partial}{\partial \eta}-c\frac{\partial}{\partial \xi}\right)\left( c\frac{\partial}{\partial \eta}-c\frac{\partial}{\partial \xi}\right)y \\
&= \left( \frac{\partial^2}{\partial \xi^2} -\frac{\partial^2}{\partial \xi\partial\eta}- \frac{\partial^2}{\partial \eta\partial\xi}  +\frac{\partial^2}{\partial \eta^2}\right)y \\
\therefore 4c^2\frac{\partial^2y}{\partial \xi\partial\eta}&=0 \\
\implies \frac{\partial y}{\partial \eta}&=G(\eta) \\
\implies y&= \int^\eta G(u)\mathrm{d}u+f(\xi)\\
&=f(x+ct)+g(x-ct) \qed
\end{align*}

Q: d'Alembert's solution
A: We are given \[
\frac{\partial ^2 y}{\partial t^2} = c^2 \frac{\partial^2 y}{\partial x^2}
\]
where $y(x,0)=\phi(x)$ and $ \frac{\partial y}{\partial t} (x,0)=\psi(x)$. Therefore, as the general solution is $y=f(x+ct)+g(x-ct)$:
\begin{align*}
\phi(x)&=f(x)+g(x) \\
\implies \phi'(x)&=f'(x)+g'(x) \tag{1}\\
\psi(x)&=cf'(x)-cg'(x) \\
\implies \frac{1}{c}\psi(x)&=f'(x)-g'(x) \tag{2} \\
\therefore f'(x)&=\frac{1}{2}\left[ \phi'(x)+\frac{1}{c}\psi(x) \right] \\
\implies f(x)&=\frac{1}{2}\phi(x)+\frac{1}{2c}\int_0^x\psi(u)\mathrm{d}u \\
\therefore g'(x)&=\frac{1}{2}\left[ \phi'(x)-\frac{1}{c}\psi(x) \right] \\
\implies g(x)&=\frac{1}{2}\phi(x)-\frac{1}{2c}\int_0^x\psi(u)\mathrm{d}u
\end{align*}
As the above holds for all $x$, we may let $x\to$ $x+ct$ and $x\to$ $x-ct$ in the two solutions respectively. This gives the d'Alembert's solution
\begin{align*}
u(x,t)&=f(x+ct)+g(x-ct) \\
&=\frac{1}{2}\phi(x+ct)+\frac{1}{2c}\int_0^{x+ct}\psi(u)\mathrm{d}u + \frac{1}{2}\phi(x-ct)-\frac{1}{2c}\int_0^{x-ct}\psi(u)\mathrm{d}u \\
&=\frac{1}{2}\left[ \phi(x+ct)+\phi(x-ct) \right]+\frac{1}{2c}\int_{x-ct}^{x+ct}\psi(u)\mathrm{d}u \\
\end{align*}
Therefore we get the particular solution without solving for $f,g$ explicitly.

Q: Deriving the diffusion equation
A: Consider a one dimensional lattice structure $x=0,\pm\Delta x\ldots$ and define the concentration $c(x,t)0$ as the \emph{expected} number of particles at $x$ at time $t$. Now assume each particle take a random walk such that it moves left with probability $p$, right with probability $p$ and stays put with probability $1-2p$. Therefore
\begin{align*}
c(x,t+\Delta t)-c(x,t)&=p [c(x+\Delta x,t)-2c(x,t)+c(x-\Delta x,t)] \\
&=p [c(x+\Delta x,t)-c(x,t)+c(x-\Delta x,t)-c(x,t)] \\
\implies \Delta t \frac{\partial c}{\partial t} +O([\Delta t]^2) &=p\left[ \Delta x \frac{\partial c}{\partial x}+ \frac{[\Delta x]^2}{2}\frac{\partial^2 c}{\partial x^2} +O([\Delta x]^3)\right] \\
&+ p\left[ -\Delta x \frac{\partial c}{\partial x}+ \frac{[\Delta x]^2}{2}\frac{\partial^2 c}{\partial x^2} +O([\Delta x]^3)\right] \\
&=p[\Delta x]^2\frac{\partial^2 c}{\partial x^2} + O([\Delta x]^3) \\
\therefore \frac{\partial c}{\partial t} +O(\Delta t) &= \frac{[\Delta x]^2}{\Delta t}\frac{\partial^2 c}{\partial x^2} +O\left(\frac{[\Delta x]^3}{\Delta t}\right)
\end{align*}
We now take a distinguished limit such that $\Delta t\to 0$, $\Delta x\to 0$ and $p\frac{[\Delta x]^2}{\Delta t}\to\kappa$ where $\kappa$ is a positive constant called the diffusion coefficient. This gives the diffusion equation \[
\frac{\partial c}{\partial t} = \kappa\frac{\partial^2 c}{\partial x^2}
\]

Q: Solving the diffusion equation
A: We will solve the diffusion equation \[
\frac{\partial \theta}{\partial t} = \kappa\frac{\partial^2 \theta}{\partial x^2}
\]
in the case of an infinitely long metal bar i.e.  $\theta(x,0)=0$ and $\theta(0,t)=H(t)$. We define a similarity variable \[
\eta =\frac{x}{2\sqrt{\kappa t}}.
\]
Applying the chain rule, we have
\begin{align*}
\frac{\partial \theta}{\partial x} &= \frac{\partial \theta}{\partial  \eta} \frac{\partial \eta}{\partial x}\\
&= \frac{1}{2\sqrt{\kappa t}} \theta'(\eta)\\
\implies\frac{\partial^2 \theta}{\partial x^2} &=\frac{1}{4\kappa t}\theta''(\eta) \tag{1}\\
\frac{\partial \theta}{\partial t} &= \frac{\partial \theta}{\partial \eta}\frac{\partial \eta}{\partial t}\\
&= -\frac{1}{2}\frac{x}{2\sqrt{\kappa}}\frac{1}{t^{3/2}} \theta'(\eta)\\
&= -\frac{\eta}{2t}\theta'(\eta) \tag{2} \\
\therefore -\frac{\eta}{2t}\theta' &= \kappa \frac{1}{4\kappa t}\theta''\\
\frac{\mathrm{d}\theta'}{\mathrm{d}t}+2\eta\theta'&=0 \\
\implies \int \frac{\mathrm{d}\theta'}{\theta'}&=-\int 2\eta \mathrm{d}\eta \\
\implies \ln \theta'&=-\eta^2+c \\
\implies \theta'&=Ae^{-\eta^2} \\
\implies \theta &= A\int_0^\eta e^{-u^2}\mathrm{d} u + \beta\\
&= \alpha\mathrm{erf}(\eta) + \beta \\
&= \alpha\mathrm{erf}\left(\frac{x}{2\sqrt{\kappa t}}\right) + \beta
\end{align*}
where $\mathrm{erf}(\eta) = \frac{2}{\sqrt{\pi}} \int_0^\eta e^{-u^2} \mathrm{d}u$ is the error function. Note that erf$(0)=0$ and $\mathrm{erf}(\eta)\to 1$ as $\eta\to \infty$. \\
Initial conditions: $\theta(0,t)=1=\beta$ so $\beta = 1$. $\theta(x,0)=0=\displaystyle\lim_{\eta\to\infty} \alpha\mathrm{erf}\left(\eta \right) +1=\alpha+1$ so $\alpha =-1$. So
\begin{align*}
\theta  &=1-\mathrm{erf}\left(\frac{x}{2\sqrt{\kappa t}}\right)\\
&= \mathrm{erfc}\left(\frac{x}{2\sqrt{\kappa t}}\right)
\end{align*}
where erfc is the complementary error function.

Q: If $a\mid bc$ and $(a, b) = 1$, then $a \mid c$.
A: \begin{proof}
  From Euclid's algorithm, there exists integers $u, v\in \mathbb{Z}$ such that $ua + vb = 1$. So multiplying by $c$, we have $uac + vbc = c$. Since $a \mid bc$, $a \mid $LHS. So $a \mid c$.
\end{proof}
\textbf{Corollary 1}  If $p$ is a prime and $p\mid ab$, then $p\mid a$ or $p\mid b$. (True for all $p, a, b$)
\begin{proof}
  Consider $(p, a) = p$ or $1$ because $p$ is a prime. If $(p, a) = p$, then $p \mid a$. Otherwise, $(p, a) = 1$ and $p \mid b$ by the theorem above.
\end{proof}
\textbf{Corollary 2} If $p$ is a prime and $p\mid n_1n_2\ldots n_k$, then $p \mid n_i$ for some $i$.
\begin{proof}
\begin{align*}
&p\mid n_1(n_2\ldots n_k) \\
\implies &p | n_1 \textrm{ or } p | (n_2\ldots n_k) \\
&\ldots \\
\implies &p | n_1 \textrm{ or } p | n_2 \textrm{ or }\ldots \textrm{ or }p | n_k \\
\end{align*}
\end{proof}

Q: Evolutionary trees
A: Let $\varphi(n)$ be the statement $n - 1$ mutants produces $n$ animals. Clearly $\varphi(1)$ is true. Given some tree with $n$ animals, remove the top mutant to get two sub-trees, with $n_1$ and $n_2$ animals, where $n_1 + n_2 = n$. If $\varphi(k)$ is true $\forall k < n$, then $\varphi(n_1)$ and $\varphi(n_2)$ are true. So the total number of mutants is $1 + (n_1 - 1) + (n_2 - 1) = n - 1$. So $\varphi(n)$ is true. Hence by strong principle of induction, $\varphi(n)$ is true for all $n$.

Q: Prove that a countable union of countable sets is countable.
A: Let $I$ be a countable index set, and for each $\alpha \in I$, let $A_\alpha$ be a countable set. $\bigcup_{\alpha\in I} A_\alpha$ is countable.
\begin{proof}
\begin{align*}
\bigcup_{\alpha\in I} A_\alpha &\leftrightarrow \bigcup_{j\geq 1} A_\alpha \textrm{ (as $I$ is countable)} \\
&\overset{\mathrm{inj}}{\to}  \mathbb{N}\times \mathbb{N} \\
&\leftrightarrow\mathbb{N}
\end{align*}
by \[
a\in \bigcup A_\alpha \mapsto (m,g_m(\alpha))
\]
where
\begin{align*}
&m = \min\{j \ | \ a\in A_{j}\}\textrm{ (exists by WOP),} \\
&g_m: A_m\to\mathbb{N} \textrm{ (because $A_m$ is countable)}.
\end{align*}
\end{proof}

Q: Techniques for solving 1st order ODEs
A: Separation of variables
Integrating factor
Exact differential
Guess
Substitution
Inversion