So I have $2$ circles with their centers at $C_1$ and $C_2$ and radii $r_1$ and $r_2$ respectively, and I have the following system of equations that solves for the $8$ points corresponding to the $4$ lines that are tangent to both circles:
$$ |\vec{P_1} - \vec{C_1}|=r_1 ~~~ (1)$$
$$ |\vec{P_2} - \vec{C_2}|=r_2 ~~~ (2)$$
$$ (\vec{P_2}-\vec{P_1})\cdot(\vec{P_1}-\vec{C_1})=0 ~~~ (3)$$
$$ (\vec{P_2}-\vec{P_1})\cdot(\vec{P_2}-\vec{C_2})=0 ~~~ (4)$$
Where $P_1$ and $P_2$ are points on the circle.
I am ultimately stuck on solving these equations.
Also if you express all the points in cartesian coordinates and solve for $(x1, y1)$ and $(x2, y2)$, then there must be $4$ different solutions for each coordinate but it gets too messy and tough to deal with, so that's why I'm asking if there's a better way of finding the solutions
Image of the situation we are dealing with: 
Here is a general and powerful method that can be applied as well to other kind of conics, e.g., ellipses.
This method uses the matrix representation of a conic section $(C)$ as follows:
$$\tag{1}ax^2+2bxy+cy^2+2dx+2ey+f=0 \ \ \iff \ \ \pmatrix{x&y&1}\underbrace{\pmatrix{a&b&d\\b&c&e\\d&e&f}}_{\Gamma}\pmatrix{x\\y\\1}=0$$
where $\Gamma$ is a so-called associated matrix to $(C)$.
The dual counterpart of $(1)$ is as follows.
A straight line with equation $ux+vy+w=0$ is tangent to $(C)$ if and only if:
$$\tag{2}Au^2+2Buv+Cv^2+2Du+2Ev+w^2=0 \ \ \iff \ \ \pmatrix{u&v&w}\underbrace{\pmatrix{A&B&D\\B&C&E\\D&E&F}}_{\Gamma^{-1}}\pmatrix{u\\v\\w}=0$$
where the associated matrix is the inverse of the previous matrix $\Gamma$.
Now, let us take particular cases, i.e., circles $C_1$ (resp $C_2$) with resp. equations $x^2+y^2-1=0$ and $x^2+y^2-8x+15=0$ (see figure below).
A straight line $ux+vy+w=0$ is thus tangent to $(C_1)$ and to $(C_2)$ if and only if its coefficients satisfy simultaneously:
$$\tag{3}\cases{u^2+v^2-w^2=0\\-15u^2+v^2-8uw-w^2=0}$$
(the computations have been done using a CAS, here Mathematica, see program below)
4 solutions are found:
$$\tag{4}\begin{cases}u=0,v=-w &\implies &y-1=0 & (L_1)\\ u=0, v=w & \implies & y+1=0 & (L_2)\\ u=w,v=\sqrt{3}w & \implies & x+\sqrt{3}y+1=0 & (L_3)\\ u=w,v=-\sqrt{3}w & \implies & x-\sqrt{3}y+1=0 & (L_4) \end{cases}$$
Remark: once you have the equations of the tangents it is easy to obtain the contact points.
Appendix: Mathematica program giving solutions $(4)$: