Let $R$ be a random vector of risky returns and let $r_f$ denote the risk free rate. Let vector of expected returns $\boldsymbol{\mu} = \operatorname{E}[R]$ and covariance matrix $\Sigma = \operatorname{Cov}(R)$.
The maximum Sharpe ratio portfolio among risky assets is called the tangency portfolio.
Quick method to tangency portfolio
Let's find the variance-frontier among ALL assets (including the risk free security) in excess return space. (The return of any zero cost portfolio, i.e. one return minus another, is an excess return.) The tangency portfolio will be the mean-variance efficient portfolio (i.e. has minimum variance given any expected return) with 0 weight on the risk free rate.
Let $x_i$ denote the weight on excess return $R_i - r_f$ (hence there's a weight of $1 - \sum_i x_i$ on the risk free rate).
\begin{equation}
\begin{array}{*2{>{\displaystyle}r}}
\mbox{minimize (over $\mathbf{x}$)} & \frac{1}{2}\mathbf{x}'\Sigma\mathbf{x} \\
\mbox{subject to} & \mathbf{x}'(\boldsymbol{\mu} - r_f) = c - r_f
\end{array}
\end{equation}
First order condition:
$$ \Sigma \mathbf{x} = \lambda \left( \boldsymbol{\mu} - r_f \right)$$
The portfolio on the frontier comprised only of risky securities ($\sum_i x_i = 1$) is:
$$ \mathbf{x}_\mathrm{tan} = \frac{\Sigma^{-1} \left( \boldsymbol{\mu} - r_f\right)}{\mathbf{1}' \Sigma^{-1} \left( \boldsymbol{\mu - r_f} \right)}$$
A full, brute force derivation
Let $R$ be a random vector of risky returns and let $r_f$ denote the risk free rate. Let vector of expected returns $\boldsymbol{\mu} = \operatorname{E}[R]$ and covariance matrix $\Sigma = \operatorname{Cov}(R)$.
A classic way to proceed is to first solve for the mean-variance frontier among risky assets, that is, given a desired expected return $c$, find the portfolio weights $\mathbf{w}$ which minimize variance.
\begin{equation}
\begin{array}{*2{>{\displaystyle}r}}
\mbox{minimize (over $w$)} & \frac{1}{2}\mathbf{w}'\Sigma \mathbf{w} \\
\mbox{subject to} & \boldsymbol{\mu}' \mathbf{w} = c \\
& \mathbf{1}'\mathbf{w} = 1
\end{array}
\end{equation}
This is a convex optimization problem because it has a convex objective subject to affine constraints. Furthermore Slater's condition is satisfied hence the first order conditions are necessary and sufficient for an optimum. Form the Lagrangian:
$$\mathcal{L} = \mathbf{w}'\Sigma\mathbf{w} - \lambda \boldsymbol{\mu}' \mathbf{w}- \gamma \mathbf{1}'\mathbf{w} $$
The first order condition with respect to $\mathbf{w}$:
$$ \frac{\partial \mathcal{L}}{\partial \mathbf{w}} = \Sigma \mathbf{w} - \lambda \boldsymbol{\mu} - \gamma \boldsymbol{1} = \mathbf{0}$$
Assuming $\Sigma$ is full rank (i.e. no risk free assets or linearly dependent assets):
$$ \mathbf{w} = \lambda \Sigma^{-1} \boldsymbol{\mu} + \gamma \Sigma^{-1} \mathbf{1} $$
Now we need to solve for multipliers $\lambda$ and $\gamma$. Using $\boldsymbol{\mu}'\mathbf{w} = c$ and $\mathbf{1}'\mathbf{w} = 1$
$$ \lambda \boldsymbol{\mu'} \Sigma^{-1} \boldsymbol{\mu} + \gamma \boldsymbol{\mu'} \Sigma^{-1} \mathbf{1} = c \quad \quad \lambda \mathbf{1'} \Sigma^{-1} \boldsymbol{\mu} + \gamma \mathbf{1'} \Sigma^{-1} \mathbf{1} = 1$$
Looks complicated but this is just a simple 2 equation, 2 variable linear system. Let:
$$s_{\mu\mu} = \boldsymbol{\mu'} \Sigma^{-1} \boldsymbol{\mu} \quad \quad s_{1\mu} = \boldsymbol{\mu'} \Sigma^{-1} \mathbf{1} \quad \quad s_{11} = \mathbf{1}' \Sigma^{-1} \mathbf{1}$$
The system and solution is:
$$ \begin{bmatrix} s_{\mu\mu} & s_{1\mu} \\ s_{1\mu} & s_{11} \end{bmatrix} \begin{bmatrix} \lambda \\ \gamma \end{bmatrix} = \begin{bmatrix} c \\ 1 \end{bmatrix} \quad \quad \begin{bmatrix} \lambda \\ \gamma \end{bmatrix} = \begin{bmatrix} s_{\mu\mu} & s_{1\mu} \\ s_{1\mu} & s_{11} \end{bmatrix}^{-1} \begin{bmatrix} c \\ 1 \end{bmatrix} $$
Hence:
$$ \mathbf{w} = \begin{bmatrix} \Sigma^{-1} \boldsymbol{\mu} & \Sigma^{-1} \mathbf{1} \end{bmatrix} \begin{bmatrix} s_{\mu\mu} & s_{1\mu} \\ s_{1\mu} & s_{11} \end{bmatrix}^{-1} \begin{bmatrix} c \\ 1 \end{bmatrix} $$
This shows that the mean-variance portfolio is a linear with varying weights on vectors $\Sigma^{-1} \boldsymbol{\mu}$ and $\Sigma^{-1} \mathbf{1} $. This in turn implies the so called two-fund separation result, that the mean-variance frontier is composed of linear combinations of portfolios $\frac{\Sigma^{-1} \boldsymbol{\mu}}{\mathbf{1}'\Sigma^{-1} \boldsymbol{\mu}}$ and $\frac{\Sigma^{-1}\mathbf{1}}{\mathbf{1}'\Sigma^{-1}\mathbf{1}}$. In portfolio-weight space, the mean-variance frontier is a line.
Finding the max sharpe ratio:
You can show that at an optimum $\mathbf{w}$ the variance is given by:
$$\begin{align*}
\mathbf{w}'\Sigma \mathbf{w} &= \lambda c + \gamma \\
&= \begin{bmatrix} c & 1 \end{bmatrix} \begin{bmatrix} s_{\mu\mu} & s_{1\mu} \\ s_{1\mu} & s_{11} \end{bmatrix}^{-1} \begin{bmatrix} c \\ 1 \end{bmatrix} \\
&= \frac{s_{11}c^2 - 2s_{1u}c +s_{uu}}{s_{11}s_{uu} - s_{1u}^2}
\end{align*} $$
Hence solving for the max Sharpe ratio on the mean-variance frontier can be written as:
$$ \begin{equation}
\begin{array}{*2{>{\displaystyle}r}}
\mbox{maximize (over $c$)} & \frac{c - r_f}{\sqrt{ \frac{s_{11}c^2 - 2s_{1u}c +s_{uu}}{s_{11}s_{uu} - s_{1u}^2}} }
\end{array}
\end{equation}
$$
Multiply out constant to make problem easier. We only need
$$ \frac{d}{dc} \left[\frac{c - r_f}{\sqrt{ s_{11}c^2 - 2s_{1u}c +s_{uu} }}\right] = 0$$
$$ \sqrt{ s_{11}c^2 - 2s_{1u}c +s_{uu} } - \frac{1}{2}(c - r_f) \frac{2 cs_{11} - 2 s_{1u}}{\sqrt{ s_{11}c^2 - 2s_{1u}c +s_{uu} }} = 0$$
$$ c = \frac{s_{uu} - r_f s_{1u}}{s_{1u} - r_fs_{11}} $$
Therefore for the tangency portfolio:
$$\begin{align*} \mathbf{w}_\mathrm{tan} &= \begin{bmatrix} \Sigma^{-1} \boldsymbol{\mu} & \Sigma^{-1} \mathbf{1} \end{bmatrix} \begin{bmatrix} s_{\mu\mu} & s_{1\mu} \\ s_{1\mu} & s_{11} \end{bmatrix}^{-1} \begin{bmatrix} \frac{s_{uu} - r_f s_{1u}}{s_{1u} - r_fs_{11}} \\ 1 \end{bmatrix} \\
&= \begin{bmatrix} \Sigma^{-1} \boldsymbol{\mu} & \Sigma^{-1} \mathbf{1} \end{bmatrix} \begin{bmatrix} \frac{1}{s_{1u} - r_f s_{11}} \\ \frac{-r_f}{s_{1u}- r_f s_{11}} \end{bmatrix}\\
&= \frac{ \Sigma^{-1} (\boldsymbol{\mu} - r_f)}{\mathbf{1}'\Sigma^{-1} \boldsymbol{\mu} - r_f \mathbf{1}'\Sigma\mathbf{1}}\\
&= \frac{ \Sigma^{-1} (\boldsymbol{\mu} - r_f)}{\mathbf{1}'\Sigma^{-1} \left( \boldsymbol{\mu} - r_f \right) }
\end{align*}$$
Note also the minimum variance portfolio is given by:
$$ \mathbf{w}_\mathrm{mvp} = \frac{\Sigma^{-1}\mathbf{1}}{\mathbf{1}'\Sigma^{-1}\mathbf{1}}$$
Notation note:
I use bold letters for vectors, capital letters for random variables and matrices, and lowercase letters for scalars. $\mathbf{1} = \begin{bmatrix} 1 \\ 1 \\ \ldots \\ 1 \end{bmatrix}$ and $\mathbf{0} = \begin{bmatrix} 0 \\ 0 \\ \ldots \\ 0\end{bmatrix}$.