Axes of Symmetry for a General Ellipse 
Given a general ellipse 
  $$Ax^2+Bxy+Cy^2+Dx+Ey+F=0$$ 
  where $B^2<4AC$, what are the major and minor axes of symmetry in the form $ax+by+c=0$?

It is possible of course to first work out the angle of rotation such that $xy,x,y$ terms disappear, in order to get an upright ellipse of the form $x^2/p^2+y^2/q^2=1$ and proceed from there. This may involve some messy trigonometric manipulations.
Could there be another approach perhaps, considering only quadratic/diophantine and linear equations?
Addendum
Here's a graphical implementation based on the answer by Ng Chung Tak. 
Addendum 2
Based on the answers by amd and by Ng Chung Tak, the equations for the axes are

$$\color{red}{\left(y-\frac {2AE-BD}{B^2-4AC}\right)=\frac {C-A\pm \sqrt{(A-C)^2+B^2}}B\left(x-\frac {2CD-BE}{B^2-4AC}\right)}$$

Note that 
$$\frac{C-A\pm \sqrt{(A-C)^2+B^2}}B\cdot \color{lightgrey}{\frac {C-A\mp\sqrt{(A-C)^2+B^2}}{C-A\mp\sqrt{(A-C)^2+B^2}}}=-\frac B{C-A\mp\sqrt{(A-C)^2+B^2}}$$
i.e. it is equal to the negative of its own reciprocal. Hence the equations for the axes can also be written as 

$$\color{red}{\left(x-\frac {2CD-BE}{B^2-4AC}\right)=-\frac {C-A\mp \sqrt{(A-C)^2+B^2}}B\left(y-\frac {2AE-BD}{B^2-4AC}\right)}$$

hence the two symmetrical anti-symmetrical forms for the axes. Here's the graphical implementation. 
 A: Let $f(x,y)=Ax^2+Bxy+Cy^2+Dx+Ey+F$. This function represents a quadratic surface, specifically an elliptic paraboloid in our case of $B^2<4AC$. The level curves of this surface form a family of concentric ellipses. The center of this family coincides with an extremum of the surface, so we can find it by solving the equation $\nabla F=0$. This yields $$\mathbf p=\left({BE-2CD\over4AC-B^2},{BD-2AE\over4AC-B^2}\right)$$ for the center of the family of ellipses. Alternatively, one can look for a translation that makes the linear terms vanish, which leads to the same equations, but with a bit more work.  
Translating the paraboloid clearly doesn’t change its orientation, so we need examine only the quadratic part $Q(x,y)=Ax^2+Bxy+Cy^2$, which is unchanged by translations, to determine the directions of the principal axes. A fairly simple way to do this is to compute the eigenvectors of the matrix $\bigl(\begin{smallmatrix}A&B/2\\B/2&C\end{smallmatrix}\bigr)$ of this quadratic form. It’s a symmetric real matrix, so its eigenspaces are orthogonal and correspond to the ellipses’ axes: the minor axis to the eigenspace of the greater of the two eigenvalues, the major axis to that of the lesser. A bit of work produces the eigenvalues $$\frac12\left(A+C\pm\sqrt{(A-C)^2+B^2}\right)$$ with corresponding eigenvectors $$\left(A-C\pm\sqrt{(A-C)^2+B^2},B\right).$$You can verify that these vectors are indeed orthogonal.  
Another way to compute these direction vectors is to find the extrema of $Q$ restricted to the unit circle. As with the eigenvalues, the major axis will lie in the direction of the local minima, the minor axis in the direction of the maxima. This can be solved via Lagrange multipliers (you will find that they are equal to the eigenvalues of $Q$) or by pulling $Q$ back by $\phi:t\mapsto(\cos t,\sin t)$, though this last method might lead to “messy” trigonometric manipulations.  
With the center of the ellipse and direction vectors of the axes in hand, it’s a simple matter to construct the equations of the axis lines.
A: The principal axes are:
$$\left( y-\frac{2AE-BD}{B^2-4AC} \right)=
\frac{C-A \color{red}{\pm} \sqrt{(A-C)^2+B^2}}{B}
\left( x-\frac{2CD-BE}{B^2-4AC} \right)$$
Also refer to another answer.
A: Idea for solution using multivariable calculus:
(1) Find the max an min of the functions $H(x,y) = x$ , $V(x,y) = y$ with the restriction $Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0$ (easy with Lagrange multipliers). This gives a "bounding box" for the ellipse. The center of the box will be the center of the ellipse.
(2) Find the max and min distances of $(x,y)$ to the center of the ellipse with the restriction $Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0$.
A: The trigonometric manipulations are not messy.
After centering, the equation is, in polar coordinates,
$$\rho^2=-\frac{F'}{A\cos^2\theta+B\cos\theta\sin\theta+C\sin^2\theta}.$$
The denominator
$$A\cos^2\theta+B\cos\theta\sin\theta+C\sin^2\theta=\frac12\left((A-C)\cos\frac\theta2+B\sin\frac\theta2\right)$$ achieves its extrema when
$$(A-C)\sin\frac\theta2=B\cos\frac\theta2,$$ which is elementary.

Update:
Read $\color{green}{2\theta}$, not $\color{red}{\dfrac\theta2}$ !
$$A\cos^2\theta+B\cos\theta\sin\theta+C\sin^2\theta=\frac12\left((A-C)\cos2\theta+B\sin2\theta+(A+C)\right)$$
and
$$(A-C)\sin2\theta=B\cos2\theta.$$
