Finite-difference Time-domain Method for 2D Wave Propagation [MATLAB]

function [] = FDTD()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%       FDTD.m
%
%    Author: M.V. Barnhart
%      Date: July 2016
% 
%   Description: Finite-difference time-domain method for evaluating
%                wave propagation in 2D.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    clear; clc; close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    Nx = 500; % x-data points
    Ny = 500; % y-data points
    Nt = 500; % Total number of time steps
    dt = 0.5; % Time step-size

    t = 0:dt:Nt; % Time array

    u = zeros(Nx, Ny); % Initial displacement field u(x, y)
    u_p = zeros(Nx, Ny);
    u_n = zeros(Nx, Ny);

    % Gaussian source variables
    t0 = 10;
    tp = 5;

    for k = 1:length(t)

        fprintf('Calculating displacement fields: %.2f %%',k/length(t)*100)

        for i = 2:Nx-1

            for j = 2:Ny-1    

                u_n(i,j) = 2*u_n(i,j) - u_p(i,j) + ...
                          0.5*( u(i+1,j) + u(i-1,j) - ...
                          4*u(i,j) + u(i,j+1) + u(i,j-1) );

            end

        end

        u_p = u;   % Set present displacement state to previous state
        u = u_n;   % Set next displacement state to current state

        % Gaussian source at corner of grid
        %u(Nx-1,Ny-1)= exp(-((k-t0)/tp)^2);  

        % Gaussian source at center of grid
        u(Nx/2, Ny/2) = exp(-((k - t0)/tp)^2);
        
        u(Nx/4, Ny/4) = exp(-((k - t0)/tp)^2);
        u(3*Nx/4, 3*Ny/4) = exp(-((k - t0)/tp)^2);
        u(3*Nx/4, Ny/4) = -exp(-((k - t0)/tp)^2);
        u(Nx/4, 3*Ny/4) = -exp(-((k - t0)/tp)^2);
        % Initial displacement on boundaries
        %u(1:Nx,1)=0.1*real(exp(-((k-t0)/tp)^2)); % Left boundary
        %u(1:Nx,Ny)=0.1*real(exp(-((k-t0)/tp)^2)); % Right boundary

        % Plot 2D displacement field and getframe()
        surf(u);
        shading interp;colormap('jet');
        set(gcf, 'Position', get(0, 'Screensize'));
        axis([0 Nx 0 Ny -1 1]);
        axis(gca,'off');
        A(k) = getframe(gcf);
        clc;

    end

    % Generate video from frames
    disp('Generating video...');
    v = VideoWriter('newfile12','Uncompressed AVI');
    open(v);
    writeVideo(v,A);
    close(v);
    disp('Video creation complete!');
    
end