Solved it. I made a mistake in my derivation.
Let me explain. I was trying to determine the solstice and equinox points of an orbit given its parameters such as eccentricity and longitude of ascending node. Using those orbital parameters, I calculated three orthogonal vectors which served as the coordinate system of the orbital plane (A.K.A the ecliptic). The I determined an equation for the position of the orbiting body at any point $\theta$ in the orbit. Using this and the body's axis of rotation, I derived an equation for the distance of the central body above the orbiting body's ecliptic. $D\left(\theta\right)=-r\left(\theta\right)\left[ \tilde{a}\cdot S^O_x \cos\left(\theta\right) + \tilde{a}\cdot S^O_z \sin\left(\theta\right) \right]$ A solstice occurs when this distance is at a maximum, and an equinox occurs when it's distance is zero. The equinox calculation is easy, $\theta=\text{arccot}\left(\frac{\tilde{a}\cdot S^O_x}{\tilde{a}\cdot S^O_z}\right)-\frac{\pi}{2}$ But the solstice is harder. Doing the usual approach, you find the derivative, and then set it to zero and solve for $\theta$. Problem is, you end up having multiple trigonometric functions that are difficult to reduce to one. You end up with this expression $\frac{\tilde{a}\cdot S^O_x}{\tilde{a}\cdot S^O_z}=\frac{\cos\theta+\epsilon}{\sin\theta}$ This is completely different from what I presented in the question, but that's why I prefaced by saying I made a mistake. This is deceptively simple to solve. You rewrite the expression such that $\cos\theta-\frac{\tilde{a}\cdot S^O_x}{\tilde{a}\cdot S^O_z}\sin\theta+\epsilon=0$ Then, using the identity$a\sin x + b\cos x=\sqrt{a^2+b^2}\sin\left(x+\text{sgn}(b)\arccos\left(\frac{a}{\sqrt{a^2+b^2}}\right)\right)$ we can reduce the two trigonometric functions in the expression to just one. $\frac{\sigma}{\tilde{a}\cdot S^O_z}\sin\left(\theta+\arccos\left[-\frac{\tilde{a}\cdot S^O_x}{\sigma}\right]\right)+\epsilon=0$where $\sigma=\sqrt{\left(\tilde{a}\cdot S^O_x\right)^2+\left(\tilde{a}\cdot S^O_z\right)^2}$. We're left with an expression that's easy to solve now. We solve for $\theta$ and get $\theta=\arcsin\left(-\frac{\epsilon(\tilde{a}\cdot S^O_z)}{\sigma}\right)-\arccos\left(-\frac{\tilde{a}\cdot S^O_x}{\sigma}\right)$
And there's the solution.