function [] = plotzeros(N)

   h = 2/N;
   [x,y] = meshgrid(-1:h:1,-1:h:1);
   phi = sqrt(x.^2+y.^2)-0.5;
   x = -1:h:1;
   y = -1:h:1;

   cnt = 0;
   for i = 1:N
      for j = 1:N
         if phi(i,j) == 0
            pt(cnt+1,1) = x(i);
            pt(cnt+1,2) = y(j);
            cnt = cnt+1;
         end
      end
   end

   for i = 1:N
      for j = 1:N+1
         if phi(i,j)*phi(i+1,j) < 0
            pt(cnt+1,1) = x(i)-(x(i+1)-x(i))*phi(i,j)/(phi(i+1,j)-phi(i,j));
            pt(cnt+1,2) = y(j);
            cnt = cnt+1;
         end
      end
   end

   for i = 1:N+1
      for j = 1:N
         if phi(i,j)*phi(i,j+1) < 0
            pt(cnt+1,1) = x(i);
            pt(cnt+1,2) = y(j)-(y(j+1)-y(j))*phi(i,j)/(phi(i,j+1)-phi(i,j));
            cnt = cnt+1;
         end
      end
   end

   cnt
   plot(pt(:,1),pt(:,2),'.')
   axis([-1 1 -1 1])
   axis square
   max(abs(sqrt(pt(:,1).^2+pt(:,2).^2)-0.5))

   dtheta = 2*pi/1000;
   for i = 1:1000
      theta = (i-1)*dtheta;
      morept(i,1) = 0.5*cos(theta);
      morept(i,2) = 0.5*sin(theta);
   end
   hold on
   plot(morept(:,1),morept(:,2),'r')
   hold off

end
