Nx = 128; xL = -1; xR = 1; h = (xR-xL)/Nx; yL = -1; yR = 1; Ny = round((yR-yL)/h);
x = linspace(xL-0.5*h,xR+0.5*h,Nx+2); y = linspace(yL-0.5*h,yR+0.5*h,Ny+2);
alpha = 3; ep = 0.013; ep2 = alpha*ep^2; dt = h^2*ep^2/(2*(h^2+2*ep^2)); Nt = 3400;
r = 0.5; phi = tanh((r-sqrt(x'.^2+y.^2))/(sqrt(2)*ep));
for it = 1:Nt
phi(1,:) = phi(2,:); phi(end,:) = phi(end-1,:);
phi(:,1) = phi(:,2); phi(:,end) = phi(:,end-1);
phi(2:end-1,2:end-1) = phi(2:end-1,2:end-1) ...
+dt*(phi(2:end-1,2:end-1).^(2*alpha-1).*(1-phi(2:end-1,2:end-1).^(2*alpha))/ep2 ...
+(phi(3:end,2:end-1)+phi(1:end-2,2:end-1)+phi(2:end-1,3:end)+phi(2:end-1,1:end-2) ...
-4*phi(2:end-1,2:end-1))/h^2);
analytic(it) = sqrt(r^2-2*dt*it);
figure(1); clf; hold on; box on; grid on
[c d] = contourf(x,y,phi,[0 0],'facecolor',[0.8500 0.3250 0.0980]);
axis image; axis([-0.5 0.5 -0.5 0.5]); drawnow;
radius(it) = mean(sqrt(sum(c(:,2:end).^2)));
end
figure(2); clf; hold on; box on;
plot([0:Nt]*dt, [r real(analytic)],'MarkerIndices',[1:200:3400])
plot([0:Nt]*dt, [r radius],'o-','MarkerIndices',[100:200:3400]);
axis([0 dt*Nt 0 0.5]); leg = legend('Analytic radius','Numerical radius');