Pseudo-shockwave on flat earth

pde.m

% originally from https://pastebin.com/96mF4PWX by flawr
% see from https://chat.stackoverflow.com/transcript/message/53851681#53851681

clc;clear;clf
% grid setup
n = 100; v = linspace(-1, 1, n);
dx = v(2)-v(1);
[x,y] = meshgrid(v,v);
dt = 0.5*dx;
% initial condition
f0 = exp(-((x-0.3).^2 + (y-0.3).^2)/0.05^2); % t = 0*dt
f1 = exp(-((x-0.3).^2 + (y-0.3).^2)/0.1^2);  % t = 1*dt
f2 = 0*x;
% wave velocity inside flat earth disk (c=0 corresponds to dirichlet BC)
c = x.^2 + y.^2 < 1;
c = conv2(c, ones(5)/5^2, 'same'); % make boundary smoother
% spatial finite-differnce laplace op
laplace = spdiags(ones(n, 1) * [1, -2, 1], -1:1, n, n); % 1d
laplace = (kron(eye(n), laplace) + kron(laplace, eye(n)))/dx^2; %2d

m = 1;alpha = 0.1; % smoothing for cmap scaling
for i = 1:1000;
    ma = max(abs(f1(:)));
    m = alpha * ma + (1-alpha) * m;
    imshow(f1(1:2:end, 1:2:end), []); caxis([-m, m]); drawnow;
    % wave equation
    % dt^2 * f - c^2 * laplace f = 0
    % finite difference approx
    %(f2 - 2*f1 + f0)/dt^2 = laplace * f1
    f2(:) = -f0(:) + 2*f1(:) + c(:).^2 .*dt^2 .* (laplace * f1(:));
    f0 = f1; f1 = f2;
end