 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.)
  
 % This solution uses fsolve.
 
 % First call to fsolve...
 options = optimset('Display','off','MaxIter',20);
 [tA,fSTar1,exitflag] = fsolve(@(t) MyF(t,E1,E2),[pi,pi],options);
 if exitflag ~= 1
     % No zero found in 20 iterations so we declare E1 and E2 to be
     % non-intersecting...
     x = [];
     y = [];
 else    
     % Found A = the first intersection point
     A = EllipsePoint(tA(1),E1);
     x(1) = A(1);
     y(1) = A(2);
     
     % Now find the second intersection point
     %         B = EllipsePoint(tB(1),E1) = EllispePoint(tB(2),E2).
     
     % Note that A and B on the same side of E1's major axis and on
     % the same side of E2's major axis.
     
     % How we set the initial guess tInit(1:2) depends upon 
     %   (a) How  tA(1) compares to pi. 
     %   (b) Whether we move "into E2" as we travel along E1 past point A
     %   (c) How  tA(2) compares to pi
     %   (d) Whether we move "into E1" as we travel along E2 past point A
     
     tol = .0000001;  % A tiny time step.
     
     % If we move a tiny bit further along E1, will we be inside E2?
     
     InsideE2 = Inside(EllipsePoint(tA(1)+tol,E1),E2);
     if tA(1)<pi
        if InsideE2
            % Conclude: 0 < tA(1) < tB(1) < pi
            tInit(1) = (pi+tA(1))/2;
        else
            % Conclude: 0 < tB(1) < tA(1) < pi
            tInit(1) = tA(1)/2;
        end
     else
         if InsideE2
            % Conclude: pi < tA(1) < tB(1) < 2pi
            tInit(1) = (tA(1)+2*pi)/2;
         else
            % Conclude: pi < tB(1) < tA(1) < 2pi
            tInit(1) = (pi+tA(1))/2;
         end
     end
     
     % If we move a tiny bit further along E2, will we be inside E1?
     InsideE1 = Inside(EllipsePoint(tA(2)+tol,E2),E1);
     if tA(2)<pi
        if InsideE1
            % Conclude: 0 < tA(2) < tB(2) < pi
            tInit(2) = (tA(2)+pi)/2;
        else
            % Conclude: 0 < tB(2) < tA(2) < pi
            tInit(2) = tA(2)/2;
        end
     else
         if InsideE1
            % Conclude: pi < tA(2) < tB(2) < 2pi
            tInit(2) = (tA(2)+2*pi)/2;
         else
            % Conclude: pi < tB(2) < tA(2) < 2pi
            tInit(2) = (tA(2)+pi)/2;
         end
     end
     options = optimset('Display','off','MaxIter',50);
     [tB,fStarB,exitflag] = fsolve(@(t) MyF(t,E1,E2),tInit,options);
     % Found the second intersection point B ...
     B = EllipsePoint(tB(1),E1);
     x(2) = B(1);
     y(2) = B(2);
 end
     
 function y = MyF(t,E1,E2)
 y = EllipsePoint(t(1),E1) - EllipsePoint(t(2),E2);
      
 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 alfa = Inside(P,E)
 % P is a point and E is an ellipse.
 % alfa is 1 if P is inside E.
 % alfa is 0 if P is not inside E.
 
 % Get the string length and foci of E...
 s = 2*E.a;
 f = sqrt(E.a^2 - E.b^2); cTau = cos(E.tau); sTau = sin(E.tau);
 F1 = [E.h+cTau*f; E.k+sTau*f];
 F2 = [E.h-cTau*f; E.k-sTau*f];
 % Can we push the string out to P?
 alfa = s > norm(P-F1) + norm(P-F2);

         