psi = 0.63;
psip = 0.27;
L = 150e-6;
Ts = 60;
k = 0.00158;
q0 = 4.62e3;
C0 =1.03;
gama = 5021.8;
Jw = 8.5e-6;
a = psi*L/Jw/Ts;
b = (1-psi-psip)*k*L/Jw;
e = q0./C0;
d = k*Ts*C0./q0;
%TIME SPAN
tt = 0:100:24000;
t = tt/Ts;
%STEP LENGTH
Nz = 51;
zz= linspace(0, L, Nz);
z = zz./L;
dz = z(2) -z(1);
%IC
IC_A = zeros(1, Nz);
IC_B = zeros(1, Nz);
IC = [IC_A, IC_B];
%Solve the pde using ode15s
[t, c] = ode15s(@(t,y)f(t, y, Nz, C0, a, b, e, d, dz, gama), t, IC);
% Extract solutions
C = c(:, 1:Nz);
g = c(:, Nz+1 : 2*Nz);
% Reinput BC
C(:, 1) = 1;
C(:, end) = (4*C(:, Nz-1)-C(:, Nz-2))./3;
%plot
figure(1)
plot(t, C(:, end), '--')
%ylim([0 1.2])
function dcqdt = f(t, c, Nz, C0, a, b, e, d, dz, gama)
dcqdt = zeros(length(c), 1);
dcdt = zeros(Nz, 1);
dqdt = zeros(Nz, 1);
%define value
C = c(1:Nz);
q = c(Nz+1 : 2*Nz);
%BC
C(1) =1;
C(end) = (4*C(Nz-1)-C(Nz-2))./3;
%interior for-loop
for i=2:Nz-1
dcdz(i) = (C(i+1) - C(i-1))./2./dz; % Using FDA
dcdt(i) = (-1./a).*(dcdz(i) + b.*(gama.*C(i)-q(i).*e));
dqdt(i) = d.*(gama.*C(i)-q(i).*e);
end
dcqdt = [dcdt; dqdt];
end
