function lindiff
%%% BENG 221 - Fall 2016 - 09/29/2016
%%% Homework 02 - problem 2 - MATLAB code
% diffusion constant
global D L
D = 0.1;
L = 1;
% domain
dx = 0.02; % step size in x dimension
dt = 0.001; % step size in t dimension
xmesh = 0:dx:L; % domain in x
tmesh = 0:dt:10; % domain in t
nx = length(xmesh); % number of points in x dimension
nt = length(tmesh); % number of points in t dimension
[xx,tt] = meshgrid(xmesh, tmesh); % mesh for calculation comfort
% Analytical solution
figure(1)
sol_anal = sin(pi()*xx/2/L).*exp(-D*(pi()/2/L)^2*tt);
surf(tmesh,xmesh,sol_anal', 'EdgeColor','none')
title('Analytical solution')
xlabel('t')
ylabel('x')
zlabel('u(x,t)')
% solution using finite differences
stepsize = D * dt / dx^2; % stepsize for numerical integration
sol_fd = zeros(nt, nx);
sol_fd(1, :) = sin(pi() * xmesh / 2 / L); % initial conditions
sol_fd(:, 1) = 0; % left boundary conditions; zero value
for t = 1:nt-1
    for x = 2:nx
        if x ~= nx
            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));
        else
            % right BC, zero flux
            sol_fd(t+1, x) = sol_fd(t, x) + stepsize * ...
            (2*sol_fd(t, x-1) - 2 * sol_fd(t, x));
        end
    end
end
figure(2)
surf(tmesh,xmesh,sol_fd','EdgeColor','none')
title('Finite differences')
xlabel('t')
ylabel('x')
zlabel('u(x,t)')

% solution using Matlab's built in "pdepe"
sol_pdepe = pdepe(0,@pdefun,@ic,@bc,xmesh,tmesh);
figure(3)
surf(tmesh,xmesh,sol_pdepe','EdgeColor','none')
title('Matlab pdepe')
xlabel('t')
ylabel('x')
zlabel('u(x,t)')

%%%%% ASSOCIATED FUNCTIONS %%%%%
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 L
u0 = sin(pi() * x / 2 / L);
function [pL, qL, pR, qR] = bc(xL, uL, xR, uR, t)
% Boundary conditions function
% In the form of general mixed BC (Robin BC):
% Left: pL + qL * f(x,t,u,du/dx) = 0
% Right: pR + qR * f(x,t,u,du/dx) = 0
pL = uL; qL = 0; % value left BC: 0
pR = 0; qR = 1; % flux right BC: 0