clear
% Generate initial conditions
N = 32; % size of mesh (NxN)
psi0 = zeros(N); % wave function will be initially zero
psidot0 = zeros(N); % initialize the time derivative
psidot0(N/2, N/2) = 1; % place the initial impulse in dpsi/dt
psidot0 = smooth2(psidot0, 2); % let it be a smooth impulse

% Set up the time axis
tfinal = 2*N;
Nt = 5*tfinal
t = linspace(0, tfinal, Nt);


% Solve the problem
t0=cputime;
psi = waves2d(psi0, psidot0, t);
calculation_time_seconds = cputime - t0

% Graphics
unix('rm /tmp/waves*.png');   % delete existing movie frames from /tmp.
minpsi = min(min(min(psi)));  % for scaling the plot.
maxpsi = max(max(max(psi)));  % for scaling the plot.

colormap(hot(256));
t0=cputime;
for i=1:Nt
   image(255*(psi(:,:,i) - minpsi)/(maxpsi - minpsi)) % psi scaled 0-255
   axis('equal')
   filename = sprintf('/tmp/waves_image%03i.png',i);
   print(filename);
end
image_time_seconds = cputime - t0

colormap(cool(256));
t0=cputime;
for i=1:Nt
   mesh(psi(:,:,i))
   axis([1,N,1,N,minpsi,maxpsi])
   caxis([minpsi,maxpsi])
   shading('faceted')
   filename = sprintf('/tmp/waves_mesh%03i.png',i);
   print(filename);
end
mesh_time_seconds = cputime - t0


