Coupled Ion exchange - Electrodialysis model using ode15s

Hello,

I have tried to model a system of coupled ion exchange - electrodialysis on MATLAB. The first screenshot shows the model equations. The second screenshot shows the rearranged equations that I used in MATLAB. For coding the system, I have discretized it in the x and z directions. I have given my data file, function file, and error message below. I realize that this has to do with a singular matrix. I am not very experienced in MATLAB and cannot find a solution to this problem. Any help would be appreciated.

Thanks in advance.

Data file:

%% Defining constants
W = 0.5;              %Width of column
L = 1;              %Length of column
a = 2;            %Surface area to volume ratio
delta = 0.1;     %Film thickness     
D = 3*10^-3;        %Film diffusion coefficient
v = 0.1;        %Velocity of fluid
epsilon = 0.4;      %Bed porosity
K = 0.2;        %Equilibrium constant between solid phase and interface
u = 0.1;            %Mobility
Itot = 4;           %Total current
F = 96485.3329;     %Faraday constant
cHstar = 2*10^-1;   %Interfacial concentration of counter ion (hydrogen)
qH = 2*10^-1;       %Solid phase concentration of counter ion (hydrogen)
cFeed = 0.0001;        %Feed concentration


%% Initializing Vectors for Independent variables
t0 = 0;             %Initial time
tf = 10;           %Final time
dt = 0.01;           %Time step
t = [0:dt:tf];      %Time vector
z = [0:0.1:L];      %Mesh generation for z direction
x = [0:0.1:W];      %Mesh generation for x direction
n = numel(z);       %Number of elements in z
m = numel(x);       %Number of elements in x


%% Initialization of Dependent variables
c0 = zeros(n,1) ;               %Initial bulk concentration vector
c0(1) = cFeed;                  %Initial bulk concentration at z = 0
q0 = zeros(n,m);                %Initial solid phase concentration matrix. Dependent on x and z directions.
phi0 = zeros(n,m);              %Initial electric potential matrix. Dependent on x and z directions.
J0 = zeros(n,m);                %Initial Solid phase flux matrix. Dependent on x and z directions.
y0 = [c0;reshape(q0,n*m,1)];    %Initial y vector appending together c0 and q0


[T,Y] = ode15s(@(t,y)EDIfun2(t,y,z,x,n,m,W,L,a,delta,D,v,epsilon,K,u,Itot,F,cHstar,qH,cFeed),t,y0);
plot(T,Y);

Function file:

function DyDt=EDIfun2(t,y,z,x,n,m,W,L,a,delta,D,v,epsilon,K,u,Itot,F,cHstar,qH,cFeed);
%% Variables being allocated zero vectors
c = zeros(n,1);         %Bulk concentration vector
q = zeros(n,m);         %Solid phase concentration matrix. Dependent on x and z directions.
phi = zeros(n,m);       %Potential matrix. Dependent on x and z directions.
J = zeros(n,m);         %Solid phase flux matrix. Dependent on x and z directions.
cstar = zeros(n,m);     %Cstar variable representing
DcDt = zeros(n,1);      %Time differential vector of bulk concentration
DqDt = zeros(n,m);      %Time differential matrix of solid phase concentration. Dependent on x and z directions.
DyDt = zeros(n+n*m,1);  %Time differential vector of concatenated c and q


%% Discretizing in x and z directions
DcDz = zeros(n,1);      %Initial vector for bulk concentration change with z direction
DqDx = zeros(n,m);      %Initial matrix for solid phase concentration change in x and z directions
DJDx = zeros(n,m);      %Initializing matrix for solid phase flux change in x and z directions
DphiDx = zeros(n,m);    %Initializing matrix for potential gradient change in x and z directions
c = y(1:n);             %Initial vector for bulk concentration allocated to y vector
q = reshape(y(n+1:n+n*m),n,m);       %Initial vector for solid phase concetration allocated to y vector

DcDz(1) = 0;            %Initial bulk concentration differential
DqDx(1,1) = 0;          %Initial solid phase concentration differential
J(1,1) = 0;             %Initial solid phase flux
DphiDx(1,1) = 0;        %Initial potential gradient 
for i = 2:n-1
    
    for j = 2:m-1
        DcDz(i) = (c(i)-c(i-1))/(z(i)-z(i-1));                  %Discretization of bulk concentration with z direction
        DqDx(i,j) = (q(i,j)-q(i,j-1))/(x(j)-x(j-1));        %Discretization of solid phase concentration with x and z direction
        DphiDx(i,j) = Itot/(F*(1-epsilon)*W*L*u*q(i,j));    %Equation for DphiDx 
        J(i,j) = u*q(i,j)*DphiDx(i,j);                      %Equation for solid phase flux
        DJDx(i,j) = (J(i,j)-J(i,j-1))/(x(j)-x(j-1));        %Discretization of solid phase flux with x and z direction
        cstar(i,j) = (q(i,j)*cHstar)/(K*qH);                %Equation for cstar
    end
end
DcDz(n) = 0;            %Final bulk concentration differential
DqDx(n,m) = 0;          %Final solid phase concentration differential
DJDx(n,m) = 0;          %Final solid phase flux
DphiDx(n,m) = 0;        %Final potential gradient


%% Differential equation with time
DcDt(1) = 1;            %Initial bulk concentration differential with time
DqDt(1,1) = 1;          %Initial solid phase concentration differential with time
for i = 2:n
    
    for j = 2:m
        DcDt(i) = -v*DcDz(i) + (((1-epsilon)/epsilon)*(DqDt(i,j) + DJDx(i,j)));
        DqDt(i,j) = a*D*(q(i,j)-(cstar(i,j))) - DJDx(i,j);
        
    end
end
DyDt = [DcDt;reshape(DqDt,n*m,1)];
end
        
    

Error:

Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN. 
> In ode15s (line 589)
  In EDIdatafile2 (line 38)
 
Warning: Failure at t=0.000000e+00.  Unable to meet integration tolerances without reducing the step size below the smallest value
allowed (7.905050e-323) at time t. 
> In ode15s (line 668)
  In EDIdatafile2 (line 38)