 function [x,y] = Intersection(E1,E2)
  % E1 is an ellipse and E2 is an ellipse that satisfy Property X.
  % x and y are column vectors that house their (possible) intersection points.
  % In particular, if E1 and E2 fail to intersect, then x and y are empty.
  % If x and y intersect in two points, then those points are given by (x(1),y(1))
  % and (x(2),y(2)). (Ignore the situation where E1 and E2 are tangent to
  % each other.)
  
 % Here is a solution that uses fminbnd...
 
 % Foci and string length for E2...
  
 s2 = 2*E2.a;
 f = sqrt(E2.a^2 - E2.b^2); cTau = cos(E2.tau); sTau = sin(E2.tau);
 u1 = E2.h+cTau*f; v1 = E2.k+sTau*f ;
 u2 = E2.h-cTau*f; v2 = E2.k-sTau*f;

 [tmin,dmin] = fminbnd(@(t) MinD(t,E1,u1,v1,u2,v2,s2), 0,2*pi);
    
 if dmin > 0
    % No intersection...
    x = []; y = [];
 else
    if tmin<pi
       tstar1 = fzero(@(t) MinD(t,E1,u1,v1,u2,v2,s2),[0 tmin]);
       P1 = EllipsePoint(tstar1,E1);
       tstar2 = fzero(@(t) MinD(t,E1,u1,v1,u2,v2,s2),[tmin,pi]);
       P2 = EllipsePoint(tstar2,E1);
    else
       tstar1 = fzero(@(t) MinD(t,E1,u1,v1,u2,v2,s2),[tmin,pi]);
       P1 = EllipsePoint(tstar1,E1);
       tstar2 = fzero(@(t) MinD(t,E1,u1,v1,u2,v2,s2),[tmin,2*pi]);
       P2 = EllipsePoint(tstar2,E1);
    end
    x = [P1(1);P2(1)]; 
    y = [P1(2);P2(2)];
 end
    
    
   function d = MinD(t,E,u1,v1,u2,v2,s2)
   % scalar t defines a point P on ellipse E
   % d is the sum of the distance from P to F1 and from P to F2
   % where F1 = (u1,v1) and F2 = (u2,v2).
   P = EllipsePoint(t,E);
   d = sqrt((P(1)-u1)^2 + (P(2)-v1)^2) + sqrt((P(1)-u2)^2 + (P(2)-v2)^2)-s2;
   
   function P = EllipsePoint(t,E)
   % scalar t defines a point P = [x;y] on ellipse E
   ctau = cos(E.tau); stau = sin(E.tau);
   x0 = E.a*cos(t); y0 = E.b*sin(t);
   x = E.h + ctau*x0 - stau*y0; y = E.k + stau*x0 + ctau*y0;
   P = [x;y];
   
   function z = E2E(t,E1,E2)
       
   z = EllipsePoint(t(1),E1) - EllipsePoint(t(2).E2);