**Location** of shortest distance between two skew lines in 3D?

Here's a solution, which I have also added to Wikipedia (https://en.wikipedia.org/wiki/Skew_lines#Nearest_Points), done completely using vectors without a need to solve equations.

Expressing the two lines as vectors:

Line 1: $\mathbf{v_1}=\mathbf{p_1}+t_1\mathbf{d_1}$

Line 2: $\mathbf{v_2}=\mathbf{p_2}+t_2\mathbf{d_2}$

The cross product of $\mathbf{d_1}$ and $\mathbf{d_2}$ is perpendicular to the lines.

$\mathbf{n}= \mathbf{d_1} \times \mathbf{d_2}$

The plane formed by the translations of Line 2 along $\mathbf{n}$ contains the point $\mathbf{p_2}$ and is perpendicular to $\mathbf{n_2}= \mathbf{d_2} \times \mathbf{n}$.

Therefore, the intersecting point of Line 1 with the above mentioned plane, which is also the point on Line 1 that is nearest to Line 2 is given by

$\mathbf{c_1}=\mathbf{p_1}+ \frac{(\mathbf{p_2}-\mathbf{p_1})\cdot\mathbf{n_2}}{\mathbf{d_1}\cdot\mathbf{n_2}} \mathbf{d_1}$

Similarly, the point on Line 2 nearest to Line 1 is given by (where $\mathbf{n_1}= \mathbf{d_1} \times \mathbf{n}$)

$\mathbf{c_2}=\mathbf{p_2}+ \frac{(\mathbf{p_1}-\mathbf{p_2})\cdot\mathbf{n_1}}{\mathbf{d_2}\cdot\mathbf{n_1}} \mathbf{d_2}$

Now, $\mathbf{c_1}$ and $\mathbf{c_2}$ form the shortest line segment joining Line 1 and Line 2.

The minimum can be found with analytical method described by @Lubin in a comment. If you want to avoid explicit differentiation, you might take a shortcut: the shortest line segment between two lines is perpendicular to both lines.

Assuming $\vec A_0$ and $\vec B_0$ are points on each line, and their respective vectors are $\vec a$ and $\vec b$, so the point of one line is $$\vec A(t)=\vec A_0+t\vec a $$ and the point of the other one is $$\vec B(s)=\vec B_0+s\vec b$$

you have the segment $\overline{AB}$ must be perpendicular to $\vec a$ and to $\vec b$, which is equivalent to zero value of respective scalar products: $$\begin{cases} (\vec A(t)-\vec B(s))\cdot \vec a = 0 \\ (\vec A(t)-\vec B(s))\cdot \vec b = 0 \end{cases}$$ which results in a system of two linear equations with unknown $t,s$.

Solve it, plug $t$ and $s$ values into $\vec A(t)$ and $\vec B(s)$ definitions and you're done.

PS.
Don't forget to consider a special case of the two lines parallel.

EDIT

$$\begin{cases} (\vec A_0+t\vec a-\vec B_0-s\vec b)\cdot \vec a = 0 \\ (\vec A_0+t\vec a-\vec B_0-s\vec b)\cdot \vec b = 0 \end{cases}$$

$$\begin{cases} t(\vec a\cdot \vec a)-s(\vec a\cdot \vec b) = (\vec B_0-\vec A_0)\cdot \vec a \\ t(\vec a\cdot \vec b)-s(\vec b\cdot \vec b) = (\vec B_0-\vec A_0)\cdot \vec b \end{cases}$$