function [phi] = curvflow(N,M)

   h = 2/N;
   for i = 1:N+1
      for j = 1:N+1
         x = -1+(i-1)*h;
         y = -1+(j-1)*h;
         phi(i,j) = sqrt(x^2+y^2)-0.8;
      end
   end
   w = -1:h:1;
%   contour(w,w,phi',-0.5:0.1:0.5,'b')
%   hold on
   contour(w,w,phi',[0 0],'r')
%   hold off
   axis([-1 1 -1 1])
   axis square
   pause
 
   dt = h^2/4;
   epsilon = 1.0e-7;
   for step = 1:M
      for i = 1:N+1
         for j = 1:N+1
            phixx = (phi(min(i+1,N+1),j)-2*phi(i,j)+phi(max(i-1,1),j))/(h*h);
            phiyy = (phi(i,min(j+1,N+1))-2*phi(i,j)+phi(i,max(j-1,1)))/(h*h);
            phix = (phi(min(i+1,N+1),j)-phi(max(i-1,1),j))/(2*h);
            phiy = (phi(i,min(j+1,N+1))-phi(i,max(j-1,1)))/(2*h);
            phixy = (phi(min(i+1,N+1),min(j+1,N+1))-phi(min(i+1,N+1),max(j-1,1))-...
                     phi(max(i-1,1),min(j+1,N+1))+phi(max(i-1,1),max(j-1,1)))/(4*h*h);
%            x = -1+(i-1)*h;
%            y = -1+(j-1)*h;
%            if sqrt((x-0.6)^2+y^2) <= 0.2
%               s = 0;
%            elseif sqrt((x-0.6)^2+y^2) <= 0.4
%               s = 5*(sqrt((x-0.6)^2+y^2)-0.2);
%            else
%               s = 1;
%            end
            s = 1;
            rhs(i,j) = s*(phixx+phiyy-(phix*(phixx*phix+phixy*phiy)+...
                                       phiy*(phixy*phix+phiyy*phiy))/...
                                      (phix^2+phiy^2+epsilon^2));
         end
      end
      phi = phi+dt*rhs;
%      contour(w,w,phi',-0.5:0.1:0.5,'b')
%      hold on
      contour(w,w,phi',[0 0],'r')
%      hold off
      axis([-1 1 -1 1])
      axis square
      pause(0.01)
   end

end
