What is Cramer's rule used for? Cramer's rule appears in introductory linear algebra courses without comments on its utility.  It is a flaw in our system of pedagogy that one learns answers to questions of this kind in courses only if one takes a course on something in which the topic is used.
On the discussion page to Wikipedia's article on Cramer's rule, we find this detailed indictment on charges of uselessness, posted in December 2009.
But in the present day, we find in the article itself the assertion that it is useful for


*

*solving problems in differential geometry;

*proving a theorem in integer programming;

*deriving the general solution to an inhomogeneous linear differential equation by the method of variation of parameters;

*(a surprise) solving small systems of linear equations.  This one is what it superficially purports to be in linear algebra texts, but then elementary row operations turn out to be what is actually used.


At some point in its history, the Wikipedia article asserted that it's used in proving the Cayley–Hamilton theorem, but that's not there now.  To me the Cayley–Hamilton theorem has always been a very memorable statement, but at this moment I can't recall anything about the proof.
What enlightening expansions on these partial answers to this question can the present company offer?
 A: One place that Cramer's rule is often useful from a theoretical point of view is that if a square matrix $M$ with entries in a commutative ring $R$ has ${\rm det}(M)$ invertible in $R,$ then the inverse $M^{-1}$ still has entries in $R$. For example, a square matrix with integer entries has an inverse with integer entries if and only if its determinant is $\pm 1.$ 
I do not think this is so obvious if we pass to elementary row operations in the field of fractions of $R$ (in the case $R$ is an integral domain). Although in any (correct) computation, an integral inverse will emerge at the end if the determinant is a unit of the ring. 
A: 
Cramer's rule is helpful when proofing "Jacobi's formula", a useful
matrix derivation identity:
$$\frac{\mathrm{d}}{\mathrm{d}t} \det A(t) = \operatorname{tr} \left (\operatorname{adj}(A(t)) \, \frac{\mathrm{d}A(t)}{\mathrm{d}t}\right )$$

For further information, have a look at Wikipedia. Jacobi's formula is used in vector analysis for example for the divergence theorem. On page 2 of this document about the divergence theorem, you can see how it is done in detail.
ADDED LATER (may 2019):
Another example for the application of Cramer's rule - or at least the formula about its inverse - appears as well in calculus of variations for proving that the determinant has divergence structure and hence one can pass on weak limits. This is rather magical as it is not expected and it is most certainly not an introductory example.
Namely, the formula

$$\mathrm{det}(A)  A^{-1} = \mathrm{cof}(A)^{T},$$

with $\mathrm{cof}(A)$ being the cofactor matrix and $T$ denoting the transpose, is the starting point of the proof of this very helpful lemma and it is stated as:

Assume $n < q < \infty$ and $u_k \rightharpoonup u$ weakly in $W^{1,q}(U;\mathbb{R}^n)$.
Then
$$\mathrm{det}(Du_k) \rightharpoonup \mathrm{det}(Du) \text{ weakly in } L^{\tfrac{q}{n}}(U).$$

This lemma and its proof can be found in Evans' book "Partial Differential Equations" in Chapter 8.2.4.
A: In control theory, the closely related rule
$$A^{-1} = \frac{1}{\det A} \operatorname{Adj}(A)$$
where $\operatorname{Adj}(A)$ is the adjugate matrix, is used for going from a state space representation to a transfer function description (doing all computations by hand).
Explicitly, if we are given a state space representation of a system 
$$\begin{align}
\dot x &= Ax + Bu \\
y &= Cx + Du
\end{align}$$
and want to find the transfer function from the input $u$ to the output $y$, we can Laplace transform the equations to:
$$\begin{align}
sX &= AX + BU \\
Y &= CX + DU
\end{align}$$
and from the first equation we get $X = (sI-A)^{-1}BU$ which inserted into the second equation yields $Y = (C(sI - A)^{-1}B +D)U$ so the transfer function $G(s)$ is
$$G(s) = C(sI-A)^{-1}B + D$$
and we use the inversion rule above to calculate $(sI-A)^{-1}$.
Of course, since we can write
$$G(s) = C\frac{\operatorname{Adj}(sI-A)}{\det (sI - A)}B + D$$
we can directly observe that the poles of the system are the same as the eigenvalues of $A$, which provides information on the system's stability.
Example
We want the transfer function $G(s)$ from $u$ to $y$ in this system:
$$\begin{align}
\dot x &= \begin{pmatrix} 1 & 3 \\ 1 & 6 \end{pmatrix} + \begin{pmatrix} 1 \\ 0 \end{pmatrix} u \\
y &= \begin{pmatrix} 3 & 1 \end{pmatrix} x
\end{align}$$
so when computing $G(s)$ we have to compute:
$$(sI-A)^{-1} = \begin{pmatrix} s - 1 & -3 \\ -1 & s - 6 \end{pmatrix} = \frac{1}{s^2-7s+3} \begin{pmatrix} s - 6 & 3 \\ 1 & s - 1 \end{pmatrix}$$
so our transfer function becomes:
$$\begin{align}
G(s) &= \begin{pmatrix} 3 & 1 \end{pmatrix}\frac{1}{s^2-7s+3} \begin{pmatrix} s - 6 & 3 \\ 1 & s - 1 \end{pmatrix} \begin{pmatrix}1 \\ 0\end{pmatrix} = \\
&= \frac{1}{s^2-7s+3}\begin{pmatrix} 3 & 1 \end{pmatrix}\begin{pmatrix} s - 6 \\ 1 \end{pmatrix} = \\
&= \frac{3s-5}{s^2-7s+3}
\end{align}$$
A: If you are asking for other applications of Cramer's rule, here is an application I learnt about recently in Harold's book Galois theory:
In Harold's book, Cramer's rule is used to prove Lagrange's theroem (not the one about groups).
Let $K$ be a field, $n$ be a positive integer. The following theorem will be formulated using the obvious action of $S_n$  on $K[x_1,\ldots,x_n]$.
Lagrange's Theorem:
Let $P,Q\in K[x_1,\ldots,x_n]$, then: $\operatorname{Stab}(P)\subseteq \operatorname{Stab}(Q)$ (Stab is the stabilizer subgroup) iff there exists symmetric polynomials $s_0,s_1,\ldots,s_a,t_0,t_1,\ldots,t_b\in K[x_1,x_2,\ldots,x_n]$ such that:
$$Q=\frac{s_aP^a+s_{a-1}P^{a-1}+\cdots+s_0}{t_bP^b+t_{b-1}P^{b-1}+\cdots+t_0}$$
A: Cramer's Rule can be used for fitting polynomials to a set of points.
Say we have $n$ ordered pairs $(x_k, y_k)$ with $0<k\le n$ with no $y_k$ the same and we want to find the polynomial $P(x)$ that fits these points. Now we want to find the coefficients of this polynomial but we do know $n$ different values it holds for $n$ different values of $x$.
$$
\left(
\begin{matrix}
x_n^{n-1}&x_n^{n-2}&\cdots&x_n&1\\
x_{n-1}^{n-1}&x_{n-1}^{n-2}&\cdots&x_{n-1}&1\\
\vdots&\vdots&\ddots&\vdots&\vdots\\
x_1^{n-1}&x_1^{n-2}&\cdots&x_1&1
\end{matrix}
\right)
\left(
\begin{matrix}
a_1 \\ a_2 \\ \vdots \\a_n
\end{matrix}
\right)
=
\left(
\begin{matrix}
y_1\\y_2\\ \vdots \\ y_n
\end{matrix}
\right)
$$
Remember that a polynomial of degree $n$ has $n-1$ infliction points. On the first matrix going along we have the same value of $x_k$ at the different powers finally finishing at $0$. When it is multiplied with the column vector the coefficients of $P$ are distributed along the columns, first column being $a_1x_k^{n-1}$ and so on. Now this is a form that can directly be used with Cramer's rule and the coefficients $a_1, ...,a_n$ make up the polynomial $P(x)=\sum^n_{k=1}a_kx^{n-k}$. 
