Here is an ellipse, $\mathrm E$, whose center occurs at $x=-1$ and $y=1$ and whose semimajor axis length is $\sqrt {2/5 \,}$. Therefore, the origin is outside of the ellipse. $$\mathrm E = \{ \mathbf x \in \mathbb R^2 : 4x^2 + 3xy + 4y^2 -x +y =0 \}$$
Here it is parameterized in polar coordinates. $$4 \rho ^2 + 3 \rho^2 \cos \phi \sin \phi + 4 \rho^2 + \rho \sin \phi - \rho \cos \phi = 0 \Rightarrow \rho = \frac {\cos \phi - \sin \phi} {4 + 3\cos \phi \sin \phi}$$ Apparently $\rho \in \mathbb R \,\forall \, \phi$. But I would expect that:
$\forall \phi \in \mathcal A \subset \left [0 , 2 \pi \right ], \rho \notin \mathbb R$ ;
$\forall \phi \in \mathcal B \subset \left [0 , 2 \pi \right ]$, the "function" $\rho \left ( \phi \right ) $ should return two real numbers;
$\forall \phi \in \{ m , n \} \subset \left [0 , 2 \pi \right ], \rho \in \mathbb R$.
Why the apparent contradiction?
As I mentioned in the comments, you overlooked a couple of sources for the apparent contradictions. First, the curve is traced twice as $\phi$ ranges over $[0,2\pi)$. For half of that range $\rho$ is negative, but it is defined for all values of $\phi$. Second, when working with polar coordinates, you are dealing with rays emanating from the origin. So, when finding the values of $\rho$ that correspond to a given line through the origin, you have to split it into two opposing rays. For each ray, there is only a single value of $\rho$.
The defining equation of $\mathrm E$ is clearly satisfied by $x=y=0$, so the ellipse passes through the origin. Its center is at $(-1/5,1/5)$ and your values for the semiaxis lengths are off, too. I’m not sure what might have gone wrong in your calculation of the center, but based on your comment it appears that you used the reciprocal square roots of the eigenvalues of the quadratic part of the equation. This is correct when the equation has the form $Ax^2+Bxy+Cy^2=1$, but in general you have to apply a scale factor to those eigenvalues: after translation to eliminate the linear terms in the general equation, you have $Ax'^2+Bx'y'+Cy'^2+F'=0$, so you have to divide the quadratic form’s eigenvalues by $-F'$ when computing semiaxis lengths. If you work through this for the general equation, you’ll find that $F'=S/\Delta$, where $S$ is the determinant of the full $3\times3$ matrix of the conic and $\Delta$ is the determinant of the quadratic part—the upper-left $2\times2$ minor of the full matrix and also the product of the quadratic form’s eigenvalues.