How can I solve the generalized Lyapunov equation for an ellipse?

Is there any method or function to solve the positive definite symmetric NxN matrix P in the following matrix equation, given any NxN matrix A and the NxN symmetric matrix Q and scalars d0, d1, d2, and d3(which define the geometry of the ellipse)? d0*P+d1*(A'*P+P*A)+d2*A'*P*A+d3*(A^2'*P+P*A^2)=-Q