function [phi] = normlevel(N)

   h = 2/N;
   x = -1:h:1;
   y = -1:h:1;
   for i = 1:N+1
      for j = 1:N+1
         phi(i,j) = sqrt((x(i)+0.7).^2+y(j).^2)-0.1;
%         phi(i,j) = sqrt(x(i).^2+y(j).^2)-0.4;
      end
   end
   contour(x,y,phi',[0 0],'b')
   axis([-1 1 -1 1])
   axis square
   pause

   vn = 1;
   thetime = 0;
   step = 1;
   while thetime < 0.6
      for i = 1:N+1
         for j = 1:N+1
            if vn >= 0
               thesign = 1;
            else
               thesign = -1;
            end

            Dminus = thesign*(phi(i,j)-phi(max(i-1,1),j))/h;
            Dplus = thesign*(phi(min(i+1,N+1),j)-phi(i,j))/h;
            if Dplus >= 0 & Dminus >= 0
               Dfinalx = Dminus;
            elseif Dplus <= 0 & Dminus <= 0
               Dfinalx = Dplus;
            elseif Dminus <= Dplus
               Dfinalx = 0;
            elseif abs(Dminus) >= abs(Dplus)
               Dfinalx = Dminus;
            else
               Dfinalx = Dplus;
            end

            Dminus = thesign*(phi(i,j)-phi(i,max(j-1,1)))/h;
            Dplus = thesign*(phi(i,min(j+1,N+1))-phi(i,j))/h;
            if Dplus >= 0 & Dminus >= 0
               Dfinaly = Dminus;
            elseif Dplus <= 0 & Dminus <= 0
               Dfinaly = Dplus;
            elseif Dminus <= Dplus
               Dfinaly = 0;
            elseif abs(Dminus) >= abs(Dplus)
               Dfinaly = Dminus;
            else
               Dfinaly = Dplus;
            end

            rhs(i,j) = -vn*sqrt(Dfinalx^2+Dfinaly^2);
         end
      end

      dt = h/(2*abs(vn));
      if thetime+dt > 0.6
         dt = 0.6-thetime;
         thetime = 0.6;
      else
         thetime = thetime+dt;
      end
      phi = phi+dt*rhs;
%      quiver(x,y,v1',v2')
%      hold on
      contour(x,y,phi',[0 0],'b')
%      hold off
      axis([-1 1 -1 1])
      axis square
      pause(0.1)
      step = step+1;
   end
   thetime

end
