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]$$
Yes, it is doable. Recall your first equation,
$$x^2+y^2+z^2 = r^2$$
More generally, define,
$$x_k^2+y_k^2+z_k^2 = s_k$$
Expand your first $i$ equation,
$$\begin{aligned}
i/i_1\;&= (x-x_1)^2+(y-y_1)^2+(z-z_1)^2\\
&=\color{brown}{(x^2+y^2+z^2)}+\color{brown}{(x_1^2+y_1^2+z_1^2)}-2x_1 x-2y_1 y-2z_1 z\\
&=\color{brown}{r^2}+\color{brown}{s_1}-2x_1 x-2y_1 y-2z_1 z
\end{aligned}$$
In short, expanding your three $i$ equations, you will get a system only linear in $x,y,z$,
$$i/i_1 = r^2+s_1-2x_1\,\color{blue}x-2y_1\,\color{blue}y-2z_1\,\color{blue}z\tag1$$
$$i/i_2 = r^2+s_2-2x_2\,\color{blue}x-2y_2\,\color{blue}y-2z_2\,\color{blue}z\tag2$$
$$i/i_3 = r^2+s_3-2x_3\,\color{blue}x-2y_3\,\color{blue}y-2z_3\,\color{blue}z\tag3$$
I'm sure you know how to solve for $x,y,z$ in those three. After doing so, substitute the expressions into,
$$x^2+y^2+z^2 = r^2\tag4$$
and you will end up with a quadratic in your last unknown $i$ expressed solely in the arbitrary values $r,\, i_k,\, x_k,\, y_k,\, z_k$.
Best Answer
A numerical way to solve this would be to use the Newton-Raphson method. This method can be extended to 3 dimensions as follows:
$$\vec{i_{n+1}}=\vec{i_n}-J^{-1}(i_n)\vec{f}(i_n)$$
Where $J$ is the Jacobian matrix of the system:
$$J= \begin{bmatrix} 3i_1^2L_1+K & Z_n & Z_n \\ Z_n & 3i_2^2L_2+K & Z_n \\ Z_n & Z_n & 3i_3^2L_3+K \\ \end{bmatrix} $$
Choose an initial "guess" $\vec{i_0}$, and repeat this process. Since it's an iterative process, the more times you evaluated it, the closer you get to the solution.