The bivariate normal density can expresed as follows
p(x)=$\frac{1}{\sqrt{2\pi}\|\Sigma|}exp(-\frac{1}{2}(x-\mu)^T\Sigma^{-1}(x-\mu))$
where $\mu$ is mean vector and $\Sigma$ is variance-covariance matrix.
First you need to understand what is $\Sigma$. This is a matric of which $i,j^{th}$ element is the covariance between $i^{th}$ and $j^{th}$ variable. As in univariate case the shape of distribution depends upon variance $\sigma$, in bivariate case shape depends upon $\Sigma$. In univariate case density at a point depends upon distance of that point from $\mu$ with respect to $\sigma$. The same thing happens in bivariate case but here you can move away from mean in infinite directions. If you move in a direction in which variance is low your density would decrease faster compared with if you move in a direction in which variance is high.
For example suppose the variance of first variable is $\sigma_1^{2}$ and variance of second variable is $\sigma_2^{2}$. And the covarince between variable 1 and variable 2 is zero. Then you would get an ellitical density depending upon the values of $\sigma_1^{2}$ and $\sigma_2^{2}$. If $\sigma_1^{2}$ > $\sigma_2^{2}$; then the major axis of ellipse would be in the direction of variable 1 and minor axis ib the direction of variable 2.
Now that you understand bivariate normal density you can write likelihood function as in single variable case.
L($\mu, \Sigma$) = $\Pi_i p(x_i)$