Packed Bed Absorber - Axial Dispersion Model

% Initialisation
close all
clear all
clc
%System Set-up %
%Define Variables
D = 0.006;    % Axial Dispersion coefficient
v = 1.061;    % Superficial velocity
epsilon = 0.3;  % Voidage fraction
k = 0.1;    % Mass Transfer Coefficient
cFeed = 2;     % Feed concentration
H=2;
L = 10;         % Column length
t0 = 0;        % Initial Time
tf = 50;      % Final time
dt = 0.006;      % Time step
dx=0.01;
z = [0:0.01:L]; % Mesh generation
t = [t0:dt:tf];% Time vector
n = numel(z);  % Size of mesh grid
%Initial Conditions / Vector Creation
c0 = zeros(n,1);     
c0(1) = cFeed;
q0 = zeros(n,1);    % t = 0, q = 0 for all z, 
y0 = [c0 ; q0];     % Appends conditions together
%ODE15S Solver
[T, Y] = ode15s(@(t,y) MyFun(t,y,z,n),t,y0);
function DyDt=MyFun(~, y, z, n)
% Defining Constants
D = 0.006;    % Axial Dispersion coefficient
v = 1.061;    % Superficial velocity
epsilon = 0.3;  % Voidage fraction
k = 0.1;    % Mass Transfer Coefficient
H=2;
dx=0.01;
% Variables being allocated zero vectors
c = zeros(n,1);
q = zeros(n,1);
DcDt = zeros(n,1);
DqDt = zeros(n,1);
DyDt = zeros(2*n,1);
DcDz = zeros(n,1);
D2cDz2 = zeros(n,1);
c = y(1:n);
q = y(n+1:2*n);
for i=2:n-1
  DcDz(i) = (c(i+1)-c(i-1))/(2*dx);
  D2cDz2(i)= (c(i-1)-(2*c(i))+c(i+1))/(dx^2) ;
end
% Calculate DcDz and D2cDz2 at z=L for boundary condition dc/dz = 0
% Standard setting for q 
DqDt(1) = k*((H*c(1))-q(1));
DcDt(1) = 0.0;
% Set time derivatives in remaining points
for i=2:n
    DqDt(i) = k*((H*c(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)-2.33*DqDt(i);
end
% Concatenate vector of time derivatives
DyDt = [DcDt;DqDt];
end

I dont know whats wrong with my code. Any help will be appreciated

BC C=2 @ t>0 @z=0

IC @ t=0 &Z= 0 to L ----> c=0 q=0

[Torsten]

You only calculated DcDz and D2cDz2 for i=2,...,n-1. In your time-derivative calculation, you need DcDz and D2cDz2 for i=n which does not exist.