function lindiff
%%% Homogeneous PDE: Linear (1-D) Diffusion
%%% BENG 221 example, 9/29/2011, updated 11/4/2017

% diffusion constant
global D
D = 0.001; 

% domain
global dx
dx = 0.02; % step size in x dimension
dt = 0.1;  % step size in t dimension
xmesh = -1:0.02:1; % domain in x
tmesh = 0:0.1:10;  % domain in t


% solution using finite differences (see Week 5 class notes)
nx = length(xmesh); % number of points in x dimension
nt = length(tmesh); % number of points in t dimension
stepsize = D * dt / dx^2; % stepsize for numerical integration
sol_fd = zeros(nt, nx);
sol_fd(1, :) = (xmesh == 0)/dx; % initial conditions; delta impulse at center
sol_fd(:, 1) = 0; % left boundary conditions; zero value
sol_fd(:, nx) = 0; % right boundary conditions; zero value
for t = 1:nt-1
    for x = 2:nx-1
        sol_fd(t+1, x) = sol_fd(t, x) + stepsize * ...
            (sol_fd(t, x-1) - 2 * sol_fd(t, x) + sol_fd(t, x+1));
    end
end

figure(1)
surf(tmesh,xmesh,sol_fd')
title('Finite differences')
xlabel('t')
ylabel('x')
zlabel('u(x,t)')


% solution using Matlab's built in "pdepe"
% See: http://www.mathworks.com/help/techdoc/ref/pdepe.html
% Also: https://mse.redwoods.edu/darnold/math55/DEProj/sp02/AbeRichards/paper.pdf
help pdepe
sol_pdepe = pdepe(0,@pdefun,@ic,@bc,xmesh,tmesh);

figure(2)
surf(tmesh,xmesh,sol_pdepe')
title('Matlab pdepe')
xlabel('t')
ylabel('x')
zlabel('u(x,t)')


% function definitions for pdepe:
% --------------------------------------------------------------

function [c, f, s] = pdefun(x, t, u, DuDx)
% PDE coefficients functions

global D
c = 1;
f = D * DuDx; % diffusion
s = 0; % homogeneous, no driving term

% --------------------------------------------------------------

function u0 = ic(x)
% Initial conditions function

global dx
u0 = (x==0)/dx; % delta impulse at center

% --------------------------------------------------------------

function [pl, ql, pr, qr] = bc(xl, ul, xr, ur, t)
% Boundary conditions function

pl = ul; % zero value left boundary condition
ql = 0;  % no flux left boundary condition
pr = ur; % zero value right boundary condition
qr = 0;  % no flux right boundary condition
