% Transient 1D diffusion, uniform layer, one side impermeable,
% as in a medical drug patch, for example.
% Demonstrate a finite-difference method of numerical solution

% ALSO SEE notes attached starting p. 5

% ReactorLab.net

% SUGGSTIONS: Vary n and lambda values to get good approximation to 
% analytical solution for this case, which has an analytical solution.

% Numerical methods are most useful for problems not having available
% analytical solution, e.g., conc at layer face varies unpredictably
% or complex variation of D or initial conc within layer.

% OTHER THINGS TO DO: compute rate of diffusion out of layer vs. time; 
% compute fraction diffusing component left in layer vs. time.
 
clear all 
fprintf('------------------- \n') % separates runs in command window

% BEGIN USER INPUTS

D = 1.0e-5; % (m2/s), diffusion coefficient
cI = 1; % (mol/m3), initial concentration
c0 = 0; % new conc at layer "face" x = 0 for t>0
L = 0.01; % (m), length of layer, impermeable on one side, "interior"
n = 20; % integer number divisions of L
tfinal = 10; % (s), final time

% lambda = D*dt/dx^2 is a dimensionless diffusion coefficient
% require lambda <= 0.25 for numerical stability

lambda = 0.05;

% END USER INPUTS

% compute time step size dt (s) to get desired lambda
dt = lambda*(L/n)^2/D % where dx = (L/n), length step

% Matlab array indexes are integers starting at 1
% c is 2D concentration array that saves all position and time data
% c(i,j), where i is position index, and j is time index
% c(1,1) is at x = 0, t = 0
% c(n+1,1) is at x = L, t = 0

% specify initial conditions

c(1:n+1,1) = cI; 
t(1) = 0;
j = 1; % time index value at t = 0

jmax = 1+round(tfinal/dt); % max value of time index j
tfinal = (jmax-1) * dt; % recompute tfinal using integer jmax
m = n+1; % index for center node to match notes

% step in time

while t(j) <= tfinal
    
    % exterior node
    c(1,j+1) = c0;
    
    % interior nodes
    for i = 2:n
        c(i,j+1) = c(i,j) + lambda*(c(i+1,j)-2*c(i,j)+c(i-1,j));
    end
    
    % center node at m = n+1
    c(m,j+1) = c(m,j) + lambda*(2*c(m-1,j)-2*c(m,j)); 
    
    t(j+1) = t(j) + dt; % update time array 
    j = j+1; % increment time index
    
end

fprintf('number of time steps = %i \n',jmax-1)

fprintf('interior conc at final time = %6.4f \n',c(m,jmax))

x = linspace(0,L,n+1); % position in layer for plot

plot(x,c(1:n+1,jmax));
hold on

plot(x,c(1:n+1,round( 0.01 * jmax)));
plot(x,c(1:n+1,round( 0.1 * jmax)));
plot(x,c(1:n+1,round( 0.5 * jmax)));

tt = sprintf('Transient 1D diffusion, conc profiles at 0.01, 0.1, 0.5, 1.0 of tfinal = %4.2f (s)',tfinal);

title(tt);
ylabel('concentration (mol/m3)')
xlabel('position in layer (m)')
hold off

figure(2)
plot(t,c(m,1:jmax))
axis([0 10 0 1])
title('Transient 1D diffusion, interior conc vs. time')
ylabel('interior conc (mol/m3)')
xlabel('time (s)')