function h = plotgaussians3(p1,mu1,sigma1,p2,mu2,sigma2,p3,mu3,sigma3, ...
    xrange,yrange,linespecs,contourlines)
  % mui,sigmai each specify a 2-dim Gaussian
  % pi are their prior probabilities
  % Plots the contours for each of the Gaussians (disregarding the priors),
  % as well as the decision boundary between them (taking into acount the priors).
  % All other arguments are optional: 
  % xrange=[fromx tox] specifies range of x values to plot (default -10 10)
  % yrange=[fromy toy] specifies range of y values to plot (default same as x)
  % linespecs is a cell array specifying line specifications for (in this order)
  %   the decision boundary, and the Gaussians
  % contourlines specifies the contour line values for the Gaussians
  if nargin<12
      linespecs = {'-k',':r','--b','-.g'};
  end
  if nargin<10
      xrange = [-10 10];
  end
  if nargin<11
      yrange = xrange;
  end
  if nargin<13
      contourlines = 5 ;
  end
  [x,y] = meshgrid(fillrange(xrange),fillrange(yrange));
  g1 = gaussian2d(x,y,mu1,sigma1);
  g2 = gaussian2d(x,y,mu2,sigma2);
  g3 = gaussian2d(x,y,mu3,sigma3);
  if length(contourlines)==1
      contourlines=0:(1.5/contourlines):1.5 ;
  end
  h = figure;
  hold;
  contour(x,y,g1*sqrt(det(sigma1)),contourlines,linespecs{2});
  contour(x,y,g2*sqrt(det(sigma2)),contourlines,linespecs{3});
  contour(x,y,g3*sqrt(det(sigma3)),contourlines,linespecs{4});
  % Code for ploting the decision boundary goes here.
  % Use contour(x,y,z,[0 0],linespecs{1}) to plot the function z==0
  % as a function of x,y. You might need to make several such plots.

function filledrange = fillrange(range)
  filledrange = range(1):((range(2)-range(1))/100):range(2);
