% Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0
 % This m file plots the ellipse above.  You will
 % have to change A,B,C,D,E,F for your particular ellipse.
 %
 % theta is the angle the ellipse will be rotated in a counterclockwise
 % direction from the x-axis. The center of the ellipse will be (xc,yc)
 % the length of one semi-axis is a and the other one is b.
 %
 % The parametric form of the ellipse is then given as matrix*vector + vector
 % below:
 %
 %       x      c      -s       a*cos(tau)        xc
 %           =              *                 +
 %       y      s       c       b*sin(tau)        yc
 %
 %
 A=1.0;
 B=2.0;
 C=4.0;
 D=3.0;
 E=6.0;
 F=-1.0;
 if B == 0
    theta=0;
 else
    theta=.5*acot((A-C)/B);
 end
 coc=cos(theta);
 sic=sin(theta);
 Ap=A*coc*coc+B*coc*sic+C*sic*sic;
 Bp=B*(coc*coc-sic*sic) + 2*(C-A)*sic*coc;
 Cp=A*sic*sic - B*sic*coc + C*coc*coc;
 Dp=D*coc+E*sic;
 Ep=-D*sic+E*coc;
 Fp=F;
 T=Ep*Ep/(4*Cp) + Dp*Dp/(4*Ap) - Fp;
 xc=-Dp*coc/(2*Ap) + Ep*sic/(2*Cp);
 yc=-Dp*sic/(2*Ap) - Ep*coc/(2*Cp);
  a=sqrt(T/Ap);
  b=sqrt(T/Cp);
%
 tau=(0:2*50*pi)'/50;
 ell(1,:)=coc*a*cos(tau)'-sic*b*sin(tau)'+xc;
 ell(2,:)=sic*a*cos(tau)'+coc*b*sin(tau)'+yc;
 plot(ell(1,:),ell(2,:),'-')
 axis('equal')