% Solve the wave equation 
% u_tt = u_xx+u_yy
% with homogeneous Neumann BC 
% Modifed by Joe Sadow, 4/27/18
% requires the chebfun package

clear; clc; close all; 
warning off;

% Number of points in each direction and dt
n= 100;
dt = 4/n^2;
t = 0;

% Chebyshev Differentiation matrix
[D,xch] = cheb(n);
D2 = D^2;

% Old Dirichlet Conditions
%D2 = D2(2:end-1,2:end-1);

% Neumann BC
BC = -D([1 n+1],[1 n+1])\D([1 n+1],2:n);

% Grid points (tensor product)
[x,y] = meshgrid(xch);
%xi = x(2:end-1,2:end-1);
%yi = y(2:end-1,2:end-1);


% Initial condition
u0 = exp(-50*((x-0.3).^2+y.^2));
u1 = exp(-50*((x+dt-0.3).^2+(y).^2));
%u1 = u0;
%mesh(x,y,u0)
zlim([-.5 .5]), caxis([-.2 .2])
% A loop that never breaks
while t<10 % run for 10 seconds
    for k =1:300
        t = t+dt;
        % Leap-frog
        u = 2*u1-u0+dt^2*(D2*u1+u1*D2');
        % Update u0 and u1
        u0 = u1;
        u1 = u;
        u1([1 n+1],:) = BC*u1(2:n,:);  %Neumann BC for (x) bdry
        u1(:,[1 n+1]) = (BC*(u1(:,2:n))')';  %Neumann BC (y) bdry
    
        % Plot solution every 100 time steps
        if mod(k,100) == 0
            contourf(x,y,u)
            title(['t = ' num2str(t)],'fontsize',18);
            zlim([-.5 .5]), caxis([-.2 .2])
            drawnow
        end 
    end 
    
    if isnan(max(max(u)))
        fprintf('Divergence! Time =%d \n', t)
        break;
    end 
end