I have following equations to get $\theta$ in terms of $a_1,...,a_4$ $$ a_1cos(\theta)^3 + a_2cos(\theta)^2sin(\theta) + a_3cos(\theta)sin(\theta)^2 + a_4sin(\theta)^3 = 0\\ a_4cos(\theta)^3 - a_3cos(\theta)^2sin(\theta) + a_2cos(\theta)sin(\theta)^2 - a_1sin(\theta)^3 = 0\\ a_2cos(\theta)^3 - a_3sin(\theta)^3 - 3a_1cos(\theta)^2sin(\theta) - 2a_2cos(\theta)sin(\theta)^2 + 2a_3cos(\theta)^2sin(\theta) + 3a_4cos(\theta)sin(\theta)^2 = 0 $$
I can simplify this by setting $y = \sin(\theta), x = \cos(\theta)$
$$ x^2 + y^2 = 1\\ a_1x^3 + a_2x^2y + a_3xy^2 + a_4y^3 = 0\\ a_4x^3 - a_3x^2y + a_2xy^2 - a_1y^3 = 0\\ a_2x^3 + (2a_3 - 3a_1)x^2y + (3a_4- 2a_2)xy^2 - a_3y^3= 0 $$
Then I can (1) solve the first and the second equations to get $y$ in terms of $a_1,...,a4,x$, $$ y^2 = 1 - x^2\\a_1x^3+a_3x(1-x^2)+(a_4(1-x^2)+a_2x^2)y = 0 \\ => y = -\frac{a_1x^3+a_3x(1-x^2)}{a_4(1-x^2)+a_2x^2} $$ (2) solve the third one to get $x$ in terms of $a_1,...,a_4$, $$ (a_4(1-x^2)+a_2x^2)(a_4x^3+a_2x(1-x^2))+\\ (a_1(1-x^2)+a_3x^2)(a_1x^3+a_3x(1-x^2))=0\\ ⇒ (a_4+(a_2-a_4)x^2)(a_2-(a_2-a_4)x^2)+ \\ (a_1+(a_3-a_1)x^2)(a_3-(a_3-a_1)x^2)=0\\ ⇒ (a_4+b_2x^2)(a_2-b_2x^2)+(a_1+b_1x^2)(a_3-b_1x^2)=0\\ ⇒ a_2a_4+a_1a_3+(b_2^2+b_1^2)x^2(1-x^2) =0\\ ⇒ x = -\sqrt{\frac{d + \sqrt{(d(4c + d))}}{2d}},\\ \sqrt{\frac{d - \sqrt{(d(4c + d))}}{2d}},\\ -\sqrt{\frac{d - \sqrt{(d(4c + d))}}{2d}}, or\\ \sqrt{\frac{d + \sqrt{(d(4c + d))}}{2d}},\\ $$ where $b_1=a_3-a_1, b_2 = a_2-a_4, c=a_2a_4+a_1a_3, d=b_2^2+b_1^2,$ and (3) check the possible solutions with the last equation to finally find $\theta = cos^{-1}(x)$.
But I wonder if this is the right way. Is there any other trick I couldn't see?
$$a_1x^3 + a_2x^2y + a_3xy^2 + a_4y^3 = 0\\ a_4x^3 - a_3x^2y + a_2xy^2 - a_1y^3 = 0\\ a_2x^3 + (2a_3 - 3a_1)x^2y + (3a_4- 2a_2)xy^2 - a_3y^3= 0$$
Using $x^2 + y^2 = 1$ reduce the powers of $y$, i.e. substitute for $y^2$
$$a_1x^3 + a_2x^2y + a_3x(1-x^2) + a_4(1-x^2)y = 0\\ a_4x^3 - a_3x^2y + a_2x(1-x^2) - a_1(1-x^2)y = 0\\ a_2x^3 + (2a_3 - 3a_1)x^2y + (3a_4- 2a_2)x(1-x^2) - a_3(1-x^2)y= 0$$
Isolate y:
$$y(a_2x^2 + a_4(1-x^2)) + a_1x^3 + a_3x(1-x^2) = 0\\ -y(a_3x^2 + a_1(1-x^2)) + a_4x^3 a_2x(1-x^2) = 0\\ y((2a_3 - 3a_1)x^2 - a_3(1-x^2)) + a_2x^3 + (3a_4- 2a_2)x(1-x^2) = 0$$
Use this trick to eliminate $y$:
$$ yc_1 + R_1 = 0, \pm yc_2 + R_2 = 0 \rightarrow c_2R_1 \mp c_1R_2 = 0$$
First and second equations:
$$(a_3x^2 + a_1(1-x^2))(a_1x^3 + a_3x(1-x^2)) + (a_2x^2 + a_4(1-x^2))(a_4x^3 a_2x(1-x^2)) = 0$$
First and third equations:
$$((2a_3 - 3a_1)x^2 - a_3(1-x^2))(a_1x^3 + a_3x(1-x^2)) - (a_2x^2 + a_4(1-x^2))(a_2x^3 + (3a_4- 2a_2)x(1-x^2)) = 0$$
Resulting in two polynomials in $x$.
Solve each and test in the other equation. Note $|x| = |\cos(\theta)| \le 1$.