How to draw a contour map of $ f(x,y)=x^2+y^2+xy$ I have used a program to see that it is an ellipse but I want to know the process of thinking to actually draw the contour map myself. 
$x^2+y^2+xy=C$ for $C=0,1,2,3,...$
I can't seem to get it into an ellipse form. What should I do?
Thanks!

I assume that you have some knowledge of linear algebra.
The quadratic form $x^2 + 2 \frac{1}{2} xy + y^2$ can also be written as a dot product between two vectors
$$x^2 + 2 \frac{1}{2} xy + y^2 = \langle \left (\begin{array}{cc}  1 & 1/2 \\ 1/2 & 1 \end{array} \right ) \cdot \left (\begin{array}{c} x \\ y \end{array} \right ) , \left (\begin{array}{c} x \\ y \end{array} \right ) \rangle= \langle M v, v \rangle$$
Let's look a bit at the matrix of this quadratic form 
$$M = \left(\begin{array}{cc}  1 & 1/2 \\ 1/2 & 1 \end{array} \right )$$
It is a real symmetric matrix so it has real eigenvalues and it has a orthonormal basis of eigenvectors. The eigenvalues are $3/2$ and $1/2$ and we also obtain unit eigenvectors:
$$\left(\begin{array}{cc}  1 & 1/2\\ 1/2 & 1 \end{array} \right ) \cdot \left (\begin{array}{c} 1/\sqrt{2}\\ 1/\sqrt{2}\end{array} \right ) = 3/2 \cdot \left(\begin{array}{c} 1/\sqrt{2}\\ 1/\sqrt{2}\end{array} \right )\\
\left(\begin{array}{cc}  1 &1/2\\1/2 & 1 \end{array} \right ) \cdot \left (\begin{array}{c} -1/\sqrt{2}\\ 1/\sqrt{2}\end{array} \right ) = 1/2 \cdot \left(\begin{array}{c} -\frac{1}{\sqrt{2}}\\ 1/\sqrt{2}\end{array} \right )
$$
Therefore we have the equality: 
$$\left(\begin{array}{cc}  1 & 1/2\\ 1/2 & 1 \end{array} \right )  \cdot \left(\begin{array}{cc}  1/\sqrt{2} & -1/\sqrt{2} \\ 1/\sqrt{2} & 1/\sqrt{2} \end{array} \right) =\left(\begin{array}{cc}  1/\sqrt{2} & -1/\sqrt{2} \\ 1/\sqrt{2} & 1/\sqrt{2} \end{array} \right) \cdot \left(\begin{array}{cc}  3/2 & 0\\ 0 & 1/2 \end{array} \right )  $$
or 
$$\left(\begin{array}{cc}  1 & 1/2\\ 1/2 & 1 \end{array} \right )  =  \left(\begin{array}{cc}  1/\sqrt{2} & -1/\sqrt{2} \\ 1/\sqrt{2} & 1/\sqrt{2} \end{array} \right) \cdot \left(\begin{array}{cc}  3/2 & 0\\ 0 & 1/2 \end{array} \right ) \cdot \left(\begin{array}{cc}  1/\sqrt{2} & -1/\sqrt{2} \\ 1/\sqrt{2} & 1/\sqrt{2} \end{array} \right)^{-1} $$
we wrote the matrix of the quadratic form $M$ as 
$$M = R\cdot D \cdot R^{-1}$$
where $D$ is the diagonal matrix $\left(\begin{array}{cc}  3/2 & 0\\ 0 & 1/2 \end{array} \right )$ formed with the eigenvalues and $R$ is the rotation matrix formed with the corresponding unit eigenvectors $R = \left(\begin{array}{cc}  1/\sqrt{2} & -1/\sqrt{2} \\ 1/\sqrt{2} & 1/\sqrt{2} \end{array} \right) = \left(\begin{array}{cc}  \cos\frac{\pi}{4} & -\sin \frac{\pi}{4} \\ \sin \frac{\pi}{4} & \cos\frac{\pi}{4} \end{array} \right) $
The equation in the vector $v$ giving the curve is 
$$\langle Mv, v \rangle = C$$
or 
$$\langle R D R^{-1} v, v \rangle = C$$
The rotation matrix $R$ preserves the dot product $\langle \cdot, \cdot \rangle$ so we have $ \langle R D R^{-1} v, v \rangle = \langle  D R^{-1} v, R^{-1} v \rangle $
Therefore the equation of the curve is 
$$\langle  D R^{-1} v, R^{-1} v \rangle  = C$$
Denote by $w = R^{-1} v$. We have $\langle  D w, w \rangle  = C$ and this equation in $w$ give an ellipse in standard form $E_0$. Our curve consists of all $v$ so that $w = R^{-1} v \in E_0$, that is $v \in R E_0$. Therefore, our curve is the rotation of the ellipse $E_0$ by the matrix $R$. Note that $R$ is a rotation of angle $\frac{\pi}{4}$ counterclockwise. 
$E_0$ is the ellipse $3/2\cdot  x^2 + 1/2\cdot y^2 = C$. Our curve $E$ is the rotation counterclockwise of $E_0$ by angle $\frac{\pi}{4}$
In general, if the matrix of the quadratic form is $M = R_{\alpha} D R_{\alpha}^{-1}$, then the curve 
$$\langle Mv, v\rangle = C$$ 
is the rotation by angle $\alpha$ of the conic
$$\langle D w, w \rangle = C$$

Obs: It's customary that for ellipses in standard form the longer axis is along the $x$-axis. Therefore, it seems preferable to use the decomposition $M = R D R^{-1}$ where the diagonal of  $D$ has elements in increasing order.  We notice that our curves are also rotations with angle $\frac{\pi}{4}$  clockwise of ellipses $1/2\cdot x^2 + 3/2\cdot y^2 =C$. 
