First, it should be noted that position angle is not just defined by "two positions". The starting point is distinct from the ending point, resulting in a difference of $180^\circ$, depending on which one is first.
The previous answers were good, I just want to offer a different perspective/derivation. One way to define the position angle is that it's the angle from North counter-clockwise to the direction in question as measured on an orthographic projection that has your starting point as the origin (assuming the astronomy standard of "North up, East left").
The algebra behind the words begins by defining the unit vector/coordinates that define the starting point:
$$\hat{n}_0 = \left[\begin{array}{c}
\cos\alpha_0 \cos\delta_0 \\
\sin\alpha_0 \cos\delta_0 \\
\sin\delta_0
\end{array}\right],$$
with your ending position taking the same form with $0\rightarrow 1$. Orthographic projection on any vector is simply the process of projecting off the component of the vector in the direction of $\hat{n}_0$. The standard formula for this is
$$\vec{v}_\perp = \vec{v} - \hat{n}_0 (\vec{v}\cdot\hat{n}_0).\tag1$$
In principle, you could apply Equation (1) numerically with $\vec{v}$ as the North pole, and then as $\hat{n}_1$, then the dot product between the normalized results will give you the position angle.
Doing it that way is probably a mistake, though. See, most positions in astronomy, $\hat{n}_0$ and $\hat{n}_1$, will be separated by a small angle, so Equation (1) will suffer badly from loss of significance. That is why it's a good idea to continue the derivation to derive a formula that doesn't have this problem.
Examining the structure of $\hat{n}_0$ is worth doing because it will simplify the algebra tremendously. It is what you get if you start with the $x$-direction unit vector, $\hat{x}$, rotate around the $y$-axis by $\delta_0$, and then rotate around the $z$-axis by $\alpha_0$. Viewed this way, finding the two vectors that are perpendicular to $\hat{n}_0$ that we need is pretty straightforward - just apply the same rotation matrices to $\hat{y}$ and $\hat{z}$. In other words, you can read them off from the columns of the rotation matrix
\begin{align}
R &= \left[\begin{array}{ccc}
\cos\alpha_0 & -\sin\alpha_0 & 0 \\
\sin\alpha_0 & \cos\alpha_0 & 0 \\
0 & 0 & 1
\end{array}\right] \left[\begin{array}{ccc}
\cos\delta_0 & 0 & -\sin\delta_0 \\
0 & 1 & 0 \\
\sin\alpha_0 & 0 & \cos\alpha_0
\end{array}\right] \\
& = \left[\begin{array}{ccc}
\cos\alpha_0 \cos\delta_0 & -\sin\alpha_0 & -\cos\alpha_0 \sin\delta_0 \\
\sin\alpha_0 \cos\delta_0 & \cos\alpha_0 & -\sin\alpha_0 \sin\delta_0 \\
\sin\delta_0 & 0 & \cos\delta_0
\end{array}\right].\tag2
\end{align}
I chose how to apply the signs to the sine functions in the two rotation matrices so that the first column would match $\hat{n}_0$.
Let's call the second column of (2) $\hat{E}'$, and the third column $\hat{N}'$. Notice that if we rotate a standard set of $x$-$y$ axes by 90 degrees counter-clockwise, then the $x$-axis corresponds to North and the $y$ to East. Thus, we can use the standard formula for the 2-dimensional component of a vector and its polar angle if we identify $\hat{n}_1\cdot\hat{N}'$ as the $x$-component and $\hat{n}_1\cdot\hat{E}'$ as the $y$. That formula is
\begin{align}
P &= \operatorname{atan2}\left(y,\,x\right)\\
& = \operatorname{atan2}\left(\sin\delta_1\cos\delta_0 - \sin\delta_0\cos\delta_1 \cos(\alpha_1-\alpha_0),\,\cos\delta_0\sin(\alpha_1-\alpha_0)\right). \tag3
\end{align}
Because $\cos\delta > 0$ for all declinations, you could make the formula look more like the standard one from textbooks. My own instinct is to avoid using $\tan\delta$ because it diverges near the poles. This formula will work for all $\alpha$ and $\delta$, as long as the two points are distinct. The only bit that remains is to fiddle with the $x$-like argument to make it behave well, numerically, when the points are near each other. To do so, use $\sin\delta_1\cos\delta_0 = \sin(\delta_1-\delta_0) + \sin\delta_0\cos\delta_1$ and $1 - \cos(\alpha_1-\alpha_0) = 2\sin^2\left(\frac{\alpha_1 - \alpha_0}{2}\right)$ to get
$$P = \operatorname{atan2}\left(\sin(\delta_1-\delta_0) + 2\sin\delta_0\cos\delta_1 \sin^2\left(\frac{\alpha_1-\alpha_0}{2}\right),\,\cos\delta_0\sin(\alpha_1-\alpha_0)\right). \tag4$$
In principle, you'd want to investigate when it's best to use (3) or (4) numerically. In practice, I suspect that (4) will do better than (3), in terms of numerical accuracy, in the vast majority of cases that astronomers care about.