I'm going to give two answers. The first I'll derive, but it uses some linear algebra, which may not be familiar to you, and converts to "xyz" coordinates to do its work. The second uses just angles, and I'll just give you the formula and a references. Ideas needed:
I'll represent an arrow from the origin $(0,0,0)$ to a point $(a, b, c)$ by writing $[a, b, c]$, and call it a "vector". I call this "the vector corresponding to the point $(a,b,c)$", just to have a name for it.
"cross product of vectors": $[x, y, z] \times [a,b, c] = [yc - zb, za - xc, xb-ya]$. If the two vectors are nonzero and not parallel, then the cross product is a vector perpendicular to both.
"dot product of vectors": $[x, y, z] \cdot [a, b, c] = ax + yb + zc$.
"length of a vector $v = [x, y, z]$: we define $\| v \| = \sqrt{x^2 + y^2 = z^2} = \sqrt{ v \cdot v }.$
If $v$ and $w$ are nonzero vectors, then the angle $\theta$ between then satisfies $ \cos \theta = \frac{v \cdot w} {\|v \| ~ \|w\|}$.
There's a function, $atan2$, in most software packages, that takes a non-origin point point $(x, y)$ of the plane and produces $atan2(x, y)$, which is the angle, measured counterclockwise from the positive $x$-axis, to the ray from $(0,0)$ through $(x, y)$. (It's a lot like "arctan", but gives correct answers for points in the 3rd and fourth quadrant by checking the sign of $x$.)
A vector $[x, y, z]$ represents a point with spherical coordinates, where the "latitude coordinate" $\theta$ is $atan2(x, y)$, and the "longitude coordinate" is $\phi = \cos^{-1} z$. This assumes that the coordinates are set up so that the $xy$ plane is the equatorial plane, and the two "poles" lie at $(0,0,1)$ and $(0,0, -1)$, and that we are measuring in units that make the radius of the sphere "1". Note that the $\phi$ coordinate of a point on the equator is 0, with the phi coordinate being $\pi/2$ (or 90 degrees) a the north pole. If you are used to measuring from $0$ (at the north pole) to $\pi$ (or 180 degrees) at the south pole, you'll need a conversion.
The inverse transformation, from $(\phi, \theta)$ coordinates to $xyz$ is $(x, y, z) = (\cos \theta \sin \phi, \sin \theta \sin \phi, \cos \phi)$.
Points $(x,y,z)$ on the plane $Ax + By + Cz = 0$ (which intersects the unit sphere in a great-circle, because the plane passes through the origin, which is the center of the sphere) can also be described by saying that each point $(x, y, z)$ corresponds to a vector $v = [x, y, z]$ with the property that $n \cdot v = 0$, where $n$ is the vector $[A, B, C]$. We call $n$ "a normal vector for the plane."
I think that's it. :)
Now let's look at your diagram. I'm going to call the center of the sphere $U = (0,0,0)$, just to have a name for it.
The great circle containing $B$ and $C$ lies in a plane that's perpendicular to the vectors $b$ and $c$ for $B$ and $C$. One such vector is $n_1 = b \times c$. The same goes for the great circle through $A$ and $B$. If we say that $a$ is the vector corresponding to $A$, then one normal to the plane of that second great circle is $n_2 = a \times b$.
The cool observation here is that the angle $\theta$ that you're looking for at $B$ is exactly the angle between these two planes...which is also exactly the angle between their normal vectors! And from item 5 above, you can compute this by
$$
\theta = \cos^{-1} \left( \frac{n_1\cdot n_2}{\|n_1\|~ \|n_2 \|}\right).
$$
So your complete "formula" for computing the angle is this:
Convert each point, $A$, $B$, $C$ in polar coordinates, to rectangular, using item 8, to get coordinate triples that we treat as vectors $a$, $b$, $c$.
Compute $n_1 = b \times c$, $n_2 = a \times b$, using item 2.
Compute $\| n_1 \| $ and $\| n_2 \|$ using item 4, and $n_1 \cdot n_2$ using item 3.
Compute $\theta = \cos^{-1} \left( \frac{n_1\cdot n_2}{\|n_1\|~ \|n_2 \|}\right)$.
That wasn't so bad, was it? :)
In general, the formula for the cosine of the angle at $B$ is
$$
\cos B = \frac{\cos b - \cos a \cos c}{\sin a \sin c}
$$
where in this case, I'm using "B" to denote the angle you've drawn at point $B$, but "$b$" to indicate the angle subtended by the points $A$ and $C$ from the origin $U$, and similarly for $a$ and $c$. (You've used these lower-case letters to label three arcs, and I'm just saying that they denote not the length of the arc in miles, but rather the measure of the arc (a number between $0$ and $180$ degrees, or between $0$ and $\pi$ radians).
Computing those subtended angles? I'd be inclined to use rectangular coordinates and dot products as I did before. :)
Also: the angle at "B" is ambiguous: if you pick the angle $\pi + B$, you'll still eventually get to point $A$. If you're writing a program, you may find that the inverse cosine messes you up sometimes because of this.
Most books on celestial navigation derive these formulas, as will any book on spherical trigonometry.
There are two laws of cosines for spherical triangles:
$$\begin{align}\cos a&=\cos b\cos c+\sin b\sin c\cos A\\
\cos A&=-\cos B\cos C+\sin B\sin C\cos a\end{align}$$
Where $A$, $B$, and $C$ are the angles of the spherical triangle and $a$, $b$, and $c$ are the opposite side. Here you know one angle and two sides of the spherical triangle with vertices at the cities and the north pole:
$$\begin{align}b=\frac{\pi}2-\phi_1,&&c=\frac{\pi}2-\phi_2,&&\text{and}\quad A=\lambda_2-\lambda_1\end{align}$$
So you can apply the first of the two laws to get
$$\cos a=\sin\phi_1\sin\phi_2+\cos\phi_1\cos\phi_2\cos(\lambda_2-\lambda_1)$$
Then the distance between the two cities is $Ra$, where $R$ is the radius of the earth and $a$ is the angle we just computed in radians. Oh, but you wanted the bearing. There is also a law of sines:
$$\frac{\sin A}{\sin a}=\frac{\sin B}{\sin b}=\frac{\sin C}{\sin c}$$
So we can now get
$$\sin C=\frac{\sin A\sin c}{\sin a}$$
And $C$ is the bearing from city $1$ to city $2$ measured east of north.
Best Answer
Define the bearing angle $\theta$ from a point $A(a_1,a_2)$ to a point $B(b_1,b_2)$ as the angle measured in the clockwise direction from the north line with $A$ as the origin to the line segment $AB$.
Then,
$$ (b_1,b_2) = (a_1 + r\sin\theta, a_2 + r\cos\theta), $$
where $r$ is the length of the line segment $AB$. It follows that $\theta$ satisfies the equation
$$ \tan\theta = \frac{b_1 - a_1}{b_2 - a_2} $$
As suggested by @rogerl we can use the $\mathrm{atan2}$ function to compute $\theta$. Let
$$ \hat{\theta} = \mathrm{atan2}(b_1 - a_1, b_2 - a_2) \in (-\pi,\pi] $$
Then the bearing angle $\theta\in[0,2\pi)$ is given by
$$ \theta = \left\{ \begin{array}{ll} \hat{\theta}, & \hat{\theta} \geq 0\\ 2\pi + \hat{\theta}, & \hat{\theta} < 0 \end{array}\right. $$
Note that the equations are given in terms of Cartesian coordinates, so it is necessary to transform to screen coordinates. I believe the formula for $\hat{\theta}$ in terms of screen coordinates $(a_1,a_2)$ and $(b_1,b_2)$ is $\hat{\theta} = \mathrm{atan2}(b_1 - a_1,a_2 - b_2)$.
You could code this function in C++ as follows.