function [phi] = flowlevel(N,finaltime)

   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).^2+y(j).^2)-0.5;
      end
   end
   contour(x,y,phi',[0 0],'b')
   axis([-1 1 -1 1])
   axis square
   pause

   for i = 1:N+1
      for j = 1:N+1
         v1(i,j) = sqrt(2)/2;
         v2(i,j) = sqrt(2)/2;
      end
   end

   thetime = 0;
   step = 1;
   while thetime < finaltime
      for i = 1:N+1
         for j = 1:N+1
            if v1(i,j) < 0
               Dphi1(i,j) = (phi(min(i+1,N+1),j)-phi(i,j))/h;
            else
               Dphi1(i,j) = (phi(i,j)-phi(max(i-1,1),j))/h;
            end
            if v2(i,j) < 0
               Dphi2(i,j) = (phi(i,min(j+1,N+1))-phi(i,j))/h;
            else
               Dphi2(i,j) = (phi(i,j)-phi(i,max(j-1,1)))/h;
            end
         end
      end
      v = sqrt(v1.^2+v2.^2);
      C = max(max(v));
      dt = 10*h/(2*C);
      if thetime+dt > finaltime
         dt = finaltime-thetime;
         thetime = finaltime;
      else
         thetime = thetime+dt;
      end
      phi = phi-dt*(v1.*Dphi1+v2.*Dphi2);
      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.01)
      step = step+1;
   end
   thetime

end
