clear all; close all;
n=40; h=2/n; x=-1:h:1; y=x;     % Domain: [-1, 1] x [-1, 1]

%%%%% % Define a level set function for the geometry
%%
%%              x^2/0.5^2 + y^2/0.4^2 = 1,   an ellipse

for i=1:n+1,
  for j=1:n+1,
     phi(i,j) = sqrt( x(i)*x(i)/0.5^2 + y(j)*y(j)/0.5^2)-1; 
  end                                    
end

z=zeros(n+1,n+1); contour(x,y,phi,z);   axis('square')  % Plot the ellipse
haja = z; curv=z;

for k=1:40;

for i=1:n+1,
  for j=1:n+1,
    [grad_l, grad_r] = eno2(n,n,i,j,phi,h,h);
    U_n = -1;
    haja(i,j) = max(U_n,0)*grad_r - min(U_n,0)*grad_l;
  end
end

%dt = 0.25*h/max(max(abs(curv)));
 dt = 0.25*h;

for i=1:n+1,
  for j=1:n+1,
     phi(i,j) = phi(i,j) + dt*haja(i,j);
  end
end

contour(x,y,phi,z,'b'); axis('square')
pause(1)

end