When I attempt to run this it produces an 'error using zeros, size inputs must be integers' for the line DcDt = zeros(n,1). The code is supposed to produce a breakthrough curve. Any help is appreciated.

cFeed = 10;     % Feed concentration
L = 0.01;         % Column length
t0 = 0;        % Initial Time
tf = 90;      % Final time                             %%%%%%%%%%% Shorter time
dt = 1;      % Time step                               %%%%%%%%%%% Smaller timestep
z = 0:0.001:L; % Mesh generation                       %%%%%%%%%%% Finer mesh
t = t0:dt:tf;% Time vector
n = numel(z);  % Size of mesh grid
%Initial Conditions / Vector Creation
c0 = zeros(n,1);% t = 0, c = 0 
c0(1)=cFeed;
q0 = zeros(n,1); % t = 0, q = 0 for all z,
q0(1) = zeros; 
y0 = [c0 ; q0];     % Appends conditions together
%ODE15S Solver
[T, Y] = ode15s(@(t,y) Myfun(t,y,z,n),[t0 tf],y0);
plot(T, Y(:,n)/cFeed,'b')
xlabel('time')
ylabel('exit conc.')

Q = trapz(T,Y(:,n));
Efficiency_ads = Q/(tf*cFeed);

function DyDt=Myfun(~,y,z,n)
% Defining Constants
D = 8.76e-4;    % Axial Dispersion coefficient
v = 0.0574;    % Superficial velocity
epsilon = 0.383;  % Voidage fraction
Kldf = 1.38e-3;    % Mass Transfer Coefficient
n = 0.5396;  % Langmuir parameter
qm = 28.75;% Saturation capacity
b = 2.052e-3;
rho = 480;
%Tc = zeros(n,1);
%q = zeros(n,1);
DcDt = zeros(n,1);
DqDt = zeros(n,1);
%DyDt = zeros(2*n,1);
zhalf = zeros(n-1,1);
DcDz = zeros(n,1);
D2cDz2 = zeros(n,1);
c = y(1:n);
q = y(n+1:2*n);
% Interior mesh points
zhalf(1:n-1)=(z(1:n-1)+z(2:n))/2;
for i=2:n-1
  DcDz(i) = ((z(i)-z(i-1))/(z(i+1)-z(i))*(c(i+1)-c(i))+(z(i+1)-z(i))/(z(i)-z(i-1))*(c(i)-c(i-1)))/(z(i+1)-z(i-1));
  D2cDz2(i) = (zhalf(i)*(c(i+1)-c(i))/(z(i+1)-z(i))-zhalf(i-1)*(c(i)-c(i-1))/(z(i)-z(i-1)))/(zhalf(i)-zhalf(i-1));
end
% Calculate DcDz and D2cDz2 at z=L for boundary condition dc/dz = 0
DcDz(n) = 0; 
D2cDz2(n) = -1.0/(z(n)-zhalf(n-1))*(c(n)-c(n-1))/(z(n)-z(n-1));
% Set time derivatives at z=0
% DcDt = 0 since c=cFeed remains constant over time and is already set as initial condition
% Standard setting for q 
DqDt(1) = Kldf*((qm*b*c^1/n(1))/(1+(b*c^1/n(1))) - q(1));
DcDt(1) = 0.0;
% Set time derivatives in remaining points
for i=2:n
    %Equation: dq/dt = k(q*-q) where q* = qs * (b*c)/(1+b*c)
    DqDt(i) = Kldf*((qm*b*c^1/n(i))/(1+(b*c^1/n(i))) - q(i));
      %Equation: dc/dt = D * d2c/dz2 - v*dc/dz - ((1-e)/(e))*dq/dt
      DcDt(i) = D*D2cDz2(i) - v*DcDz(i) - ((1-epsilon)/(epsilon))*rho*DqDt(i);
  
end
% Concatenate vector of time derivatives
DyDt = [DcDt;DqDt];
end 

Answer 1

[Rik]

L is not defined. If it is smaller than 0, this will lead to z being empty, meaning that n is 0.

L=-1;
z=0:0.001:L;
numel(z)

ans = 0

After your edit, It is now apparent you have a name collision: the fourth input is called n, but so is another parameter. From the further use it seems that n is also supposed to be a vector at some point. So it is a postive integer scalar and a scalar double ('Langmuir parameter') and a vector.

Your comments don't explain enough for me to fix the bugs.

[Andre Ighomereho]

Apologies, L must have cut out when I posted the code. L is 0.01 so bigger than 0.

[Andre Ighomereho]

Ah yes that makes sense, don't know why it took so long for me to realise. Cheers!