function main()
 
   % the function whose zero level set and inner and outer approximations will be drawn
   f = inline('60-real(z).^2-1.2*imag(z).^2-0.006*(real(z)-6).^4-0.01*(imag(z)-5).^4', 'z');
 
   M=10; i=sqrt(-1); lw=2.5;
   figure(1); clf; hold on; axis equal; axis off;
 
   if  1==0
      for p=-M:M
         for q=-M:M
            z=p+i*q;
            if f(z)>0
               plot(real(z), imag(z), 'r.')
            else
               plot(real(z), imag(z), 'b.')
            end
         end
      end
   end
 
% draw the zero level set of f
   h=0.1;
   XX = -M:h:M; YY = -M:h:M;
   [X, Y] = meshgrid (XX, YY); Z = f(X+i*Y);
   [C, H] = contour(X, Y, Z, [0, 0]);
   set(H, 'linewidth', lw, 'EdgeColor', [0;0;156]/256);
 
% plot the outer polygonal curve
   Start=5+6*i; Dir=-i; Sign=-1; 
   plot_poly (Start, Dir, Sign, f, lw, [139;10;80]/256);
 
% plot the inner polygonal curve
   Sign=1; Start=4+5*i; 
   plot_poly (Start, Dir, Sign, f, lw, [0;100;0]/256);
 
%  a dummy plot to avoid a matlab bug causing some lines to appear too thin
   plot(8.5, 7.5, '*', 'color', 0.99*[1, 1, 1]);
   plot(-4.5, -5, '*', 'color', 0.99*[1, 1, 1]);
 
   saveas(gcf, 'jordan_illustration.eps', 'psc2');
 
function plot_poly (Start, Dir, Sign, f, lw, color)
 
   Current_point = Start;
   Current_dir   = Dir;
 
   Ball_rad = 0.03;
 
   for k=1:100
 
      Next_dir=-Current_dir;
 
      % from the current point, search to the left, down, and right and see where to go next
      for l=1:3
         Next_dir = Next_dir*(Sign*i);
 
         if Sign*f(Current_point+Next_dir)>=0 & Sign*f(Current_point+(Sign*i)*Next_dir) < 0
            break;
         end
      end
 
      Next_point = Current_point+Next_dir;
 
      plot([real(Current_point), real(Next_point)], [imag(Current_point), imag(Next_point)], 'linewidth', lw, 'color', color);
 
      round_ball(Current_point, Ball_rad, color'); % just for beauty, to round off some rough corners
 
      Current_dir=Next_dir;
      Current_point = Next_point;
 
   end
 
 
function round_ball(z, r, color)
   x=real(z); y=imag(z);
   Theta = 0:0.1:2*pi;
   X = r*cos(Theta)+x;
   Y = r*sin(Theta)+y;
   Handle = fill(X, Y, color);
   set(Handle, 'EdgeColor', color);