function reini

set(gcf,'rend','z');
set(gcf,'doublebuffer','on');

fF=inline('1+x+y','x','y');
dinit=inline('dcircle(p,-.5,-.5,.2)','p');
lssim(fF,dinit,5e-3,200,4,'F=1+x+y');

function lssim(fF,dinit,dt,niter,pltfreq,tit)

h=0.025;
[xx,yy]=meshgrid(-1.1:h:1.1,-1.1:h:1.1);
p=[xx(:),yy(:)];

zz=reshape(dinit(p),size(xx));

plt(xx,yy,zz,tit)
pause

for iter=1:niter
  F=fF(xx,yy);
  zz=lsstep(zz,F,h,dt);
  zz=h*reinit2dc(zz/h,0,[.2,1.5]);

  if mod(iter,pltfreq)==0, plt(xx,yy,zz,tit); end
end
pause
for iter=1:200
  zz=h*reinit2dc(zz/h,1,[.2,1.5]);
  plt(xx,yy,zz,tit);
end

function plt(xx,yy,zz,tit)

%surf(xx,yy,0*xx,zz)
clf
hold on,contour(xx,yy,zz,-2:.1:2,'c');hold off
hold on,[foo,handle]=contour(xx,yy,zz,[0,0],'k');set(handle,'linewidth',2);hold off
view(2),axis equal,axis([-1,1,-1,1])
title(tit,'fontsize',24,'fontname','times','fontangle','italic');
drawnow

function zz=lsstep(zz,F,h,dt)

[m,n]=size(zz);
Dnx=zeros(m,n); 
Dpx=zeros(m,n); 
Dny=zeros(m,n); 
Dpy=zeros(m,n); 

Dx=(zz(2:m,:)-zz(1:m-1,:))/h;
Dy=(zz(:,2:n)-zz(:,1:n-1))/h;

Dnx(2:m,:)=Dx;
Dpx(1:m-1,:)=Dx;
Dny(:,2:n)=Dy;
Dpy(:,1:n-1)=Dy;
  
delp=sqrt(max(Dnx,0).^2+min(Dpx,0).^2+ ...
          max(Dny,0).^2+min(Dpy,0).^2);
deln=sqrt(max(Dpx,0).^2+min(Dnx,0).^2+ ...
          max(Dpy,0).^2+min(Dny,0).^2);

upd=-dt*(max(F,0).*delp+min(F,0).*deln);

zz(:)=zz(:)+upd(:);