Finite Difference Explicit Method for Fick's 2nd Law

Finite Difference Explicit Method (iteration) to solve for a PDE (Fick's 2nd Law of diffusion)

%% Finite Difference Explicit Method (Iteration) with D = diffusivity: Fick's 2nd Law of Diffusion
% by Prof. Roche C. de Guzman
clear; clc; close('all');
%% Given
xi = 0; xf = 0.6; dx = 0.04; % x range and step size = dx [m]
xL = 0; xU = 0.1; % initial value x lower and upper limits [m]
ti = 0; tf = 0.05; dt = 4e-4; % t range and step size = dt [s]
ci = 2; % initial concentration value [ng/L]
cLU = 8; % initial concentration value within x lower and upper limits [ng/L]
D = 1.5; % diffusivity or diffusion coefficient [m^2/s]
%% Calculations
% Independent variables: x and t
X = xi:dx:xf; nx = numel(X); T = ti:dt:tf; nt = numel(T); % x and t vectors and their number of elements
[x,t] = meshgrid(X,T); x = x'; t = t'; % x and t matrices
% Dependent variable: c
c = ones(nx,nt)*ci; % temporary c(x,t) matrix with rows: c(x) and columns: c(t)
% Initial values and Dirichlet boundary
I = find((X>=xL)&(X<=xU)); % index of lower and upper limits
c(I,1) = cLU; % c at t = 0 for lower and upper limits
% Iteration using Taylor's Finite Difference
for j = 1:nt-1 % t counter
for i = 1:nx-2 % x counter
c(i+1,j+1) = c(i+1,j)+((dt/(dx)^2)*D*(c(i+2,j)-2*c(i+1,j)+c(i,j))); % c(x+dx,t+dt)
% Neumann boundary (zero flux): cx(xi,t)=0 and cx(xf,t)=0 -> c(xi,t)=c(xi+dx,t) and c(xf,t)=c(xf-dx,t)
c(1,j+1) = c(2,j+1); % change c(xi) = c(xi+dx)
c(nx,j+1) = c(nx-1,j+1); % change c(xf) = c(xf-dx)
end
end