A hint in your formula's favor is that if we permute the variables, $\alpha,\beta,\gamma$, you get the same formula. Thus, it is necessary that if the coefficient of $\sin\alpha \sin\beta\cos \gamma$ is $-1$, then the other two terms with two $\sin$s in them would also have to have coefficients $-1$.
A general rule is, if you are worried about signs, check with specific examples.
What happens when $\beta=-\alpha$? Then your formula would be:
$$\cos\alpha\cos(-\alpha)\cos\gamma -\sin\alpha\sin(-\alpha)\cos\gamma-\sin\alpha\cos(-\alpha)\sin\gamma - \cos\alpha\sin(-\alpha)\sin\gamma =\\
\cos\gamma(\cos^2\alpha + \sin^2\alpha) - \sin\alpha\cos\alpha\sin\gamma + \cos\alpha\sin\alpha\sin\gamma\\ = \cos \gamma = \cos(\alpha+\beta+\gamma)$$
what happens when $\alpha=\beta=45^\circ$? then $\cos(\alpha+\beta+\gamma)=-\sin\gamma$.
Since $\sin\alpha=\sin\beta=\cos\alpha=\cos\beta=\sqrt{\frac{1}{2}}$, you can plug in and check again:
$$\frac{1}{2}\cos\gamma - \frac{1}{2}\cos\gamma -\frac{1}{2}\sin\gamma -\frac{1}{2}\sin\gamma = -\sin\gamma$$
When you get to complex numbers, you'll see an elegant way to "see" this is the right formula. The formula for $\cos(\alpha+\beta+\gamma+\dots)$ will in general have terms with an even number of $\sin$ expressions, and the coefficient for each term will be $1$ if there is a multiple of $4$ $\sin$ terms, and $-1$ if there are not.
Rearrange the two equations as follows:-
$$l_2\cos(\gamma)=l_1\sin(\alpha)-l_3\sin(\beta)$$
$$l_2\sin(\gamma)=l_3\cos(\beta)-l_1\cos(\alpha)+l_4$$
Square both sides and add the two equations to eliminate the $\gamma$ term (using the identity that $\sin^2\theta+\cos^2\theta=1$) :-
$$l_2^2(\cos^2(\gamma)+\sin^2(\gamma))=l_2^2\\=l_1^2+l_3^2-2l_1l_3(\cos(\alpha)\cos(\beta)+\sin(\alpha)\sin(\beta))+2l_3l_4\cos(\beta)-2l_1l_4\cos(\alpha)$$
To solve for $\beta$, let $u(\beta)=\cos(\beta)$, so that $\sqrt{1-u(\beta)^2}=\sin(\beta)$
Thus we have
$$u(\beta)(2l_3l_4-2l_1l_3\cos(\alpha))+l_1^2+l_3^2-l_2^2-2l_1l_4\cos(\alpha)=2l_1l_3\sin(\alpha)\sqrt{1-u(\beta)^2}$$
Let us denote the following
$$f(\alpha)=2l_3l_4-2l_1l_3\cos(\alpha)$$
$$g(\alpha)=l_1^2+l_3^2-l_2^2-2l_1l_4\cos(\alpha)$$
$$h(\alpha)=2l_1l_3\sin(\alpha)$$
We end up with a quadratic in $u(\beta)$
$$u(\beta)f(\alpha)+g(\alpha)=h(\alpha)\sqrt{1-u(\beta)^2}\\\Rightarrow u(\beta)^2(f(\alpha)^2+h(\alpha)^2)+u(\beta)(2f(\alpha)g(\alpha))+g(\alpha)^2-h(\alpha)^2=0\\\Rightarrow u(\beta)=\frac{-2f(\alpha)g(\alpha)\pm\sqrt{f(\alpha)^2g(\alpha)^2-(f(\alpha)^2+h(\alpha)^2)(g(\alpha)^2-h(\alpha)^2)}}{f(\alpha)^2+h(\alpha)^2}$$
Culminating in $\beta$ being
$$\beta=\arccos\left[\frac{-2f(\alpha)g(\alpha)\pm\sqrt{f(\alpha)^2g(\alpha)^2-(f(\alpha)^2+h(\alpha)^2)(g(\alpha)^2-h(\alpha)^2)}}{f(\alpha)^2+h(\alpha)^2}\right]$$
Best Answer
Its likely that the fsolve isn't working in matlab because you're initial guess isn't close enough, causing the optimization software to get stuck. If you try starting with some values that approximate the answer you want, you might get an exact answer from matlab.
Some values that worked for me were
$A = 4000$
$\alpha = 13^\circ$
$\beta = -3^\circ$
$\gamma = -3^\circ$
hope this helps.