%Christopher Lum
%christopher.w.lum@gmail.com
%
%Investigate numerical solutions to partial differential equations (PDEs)
%
%This file is designed to accompany the YouTube video at
%https://youtu.be/ZSNl5crAvsw

%Version History
%11/06/19: Created from lecture notes to accompany YouTube
%11/10/19: Added YouTube URL

clear
clc
close all

%% Define constants
Nt = 1000;  %number of temporal steps
Nx = 100;   %number of spatial steps

tMin = 0.1;
tMax = 10;

xMin = -5;
xMax = 5;

t = linspace(tMin, tMax, Nt);
x = linspace(xMin, xMax, Nx);

%Compute secondary parameters
DeltaT = t(2) - t(1);
DeltaX = x(2) - x(1);

%% Numerically solve the PDE
%Initialize a container to hold all the p values
p = zeros(Nt, Nx);

%Set initial conditions
p(1,:) = AnalyticalSolution(t(1), x);

%plot the initial condition
figure
plot(x, p(1,:))

%iterate through time values
for i=2:Nt
    
    %iterate through spatial values
    for j=1:Nx
        %deal with boundary conditions
        if(j==1 || j==Nx)
            %at the first or last column at the edges of the boundary
            p(i,j) = 0;
        else
            %inside the middle of the matrix
            p_im1_jp1   = p(i-1,j+1);   %previous time point and forward spatial point
            p_im1_j     = p(i-1, j);    %previous time point and current spatial point
            p_im1_jm1   = p(i-1, j-1);  %previous time point and previous spatial point
            
            p(i,j) = 0.5*(DeltaT/(DeltaX^2))*(p_im1_jp1 - 2*p_im1_j + p_im1_jm1) + p_im1_j;
        end
    end
end

%% Plot the results
figure
for i=1:20:Nt
    clf
    hold on
    
    pAnalytical = AnalyticalSolution(t(i), x);
    pNumerical  = p(i, :);
    
    plot(x, pNumerical, x, pAnalytical, 'LineWidth', 2)
    xlabel('x')
    ylabel('p(t,x)')
    title(['p(t,x) @ t = ',num2str(t(i))])
    grid on
    legend('Numerical','Analytical')
    axis([min(x) max(x) 0 max(max(p))])
    drawnow
end
