Probably no longer in need of a response but...
The way I solved this was to find the "levels" that the contours need to be drawn at to contain a certain percentile of the 2D gaussian. Given "alpha" (%), this depends on critical values of the chi-square distribution, for which you can use "chi2inv". More info here.
Below an example of how to calculate the "levels" for a couple alphas. It's done by using a purely diagonal covariance matrix ("unrotated"). After obtaining "pdlevels" you can feed those to contour...
mu = [7, 3];
cov = [1 0.2; 0.2, 1];
alpha = [0.2 0.7];
nAlpha = length(alpha);
% walk out from the mean to points such that all points closer
% comprise the central alpha% of the distribution; choose the
% direction along which this distance is shortest so that we
% can return the radius
pdpoints = repmat(mu, 2*nAlpha,1) + ...
repmat(sqrt(chi2inv(alpha, 2))*...
sqrt(cov),2,1).*sortrows(repmat(eye(2), nAlpha,1));
% calculate the value of the (unrotated) gaussian there.
pdlevels = mvnpdf(pdpoints(1:nAlpha,:), mu, diag(diag(cov)));
