You can use Eigen decomposition to come up with the parameters for the ellipse. You then scale the size of the ellipse based on how many standard deviations you want to include. The components of the ellipse are computed from the eigen decomposition of your data matrix. I pretty much copied the code below from here.
%for demonstration, I generated some data: 100 points in 2D
X = normrnd(0,2,[100 2]);
%the mean should be at zero
Mu = mean(X);
X0 = bsxfun(@minus, X, Mu);
%scale the size of the ellipse based on standard deviatioon
STD = 2;
%covers around 95% of population
conf = 2*normcdf(STD)-1;
%inverse chi-squared with dof=dimensions
scale = chi2inv(conf,2);
Cov = cov(X0) * scale;
% eigen decomposition [sorted by eigen values]
[V, D] = eig( Cov );
[D, order] = sort(diag(D), 'descend');
D = diag(D);
V = V(:, order);
t = linspace(0,2*pi,100);
% unit circle
e = [cos(t) ; sin(t)];
% scale eigenvectors
VV = V*sqrt(D);
% project circle back to orig space
e = bsxfun(@plus, VV*e, Mu');
% plot cov and major/minor axes
figure(1)
plot(e(1,:), e(2,:), 'Color','b');
hold on
scatter(X(:,1),X(:,2))
hold off
