I'm reading a paper that uses a system of Blashcke products to interpolate some arbitrary data, and I want to solve this system explicitly. I have reduced the problem to the following. Let $p_1,p_2,q_1,q_2$ be fixed complex numbers in the unit disk. We may assume $p_1$ and $p_2$ are real if need be. I want to solve the system of equations
$\begin{align*} &(p_1-c^2)(\overline{d}q_1-d)=(q_1-d^2)(\overline{c}p_1-c)\\ &(p_2-c^2)(\overline{d}q_2-d)=(q_2-d^2)(\overline{c}p_2-c) \end{align*}.$
Is there a way to solve for $c$ and $d$ where both lies in the unit disk? The issue seems to be the conjugates. Blowing up this set of equations into 4 real equations seems to be incredibly messy, and converting back to a solution with $p_1,p_2,q_1$ and $q_2$ seems to be near impossible. Polar coordinates doesn't seem to be helpful with the addition in the equations. Also, it would be helpful if could solve the system when either $c$ or $d$ real (but not both). I would appreciate a strategy for solving this type of system, not a solution. I have access to computer algebra systems if need be.