Hi everyone, I am trying to solve a Diffusion-​Convection​-Reaction PDE (1D) using Crank-Nicholson Method given by dC/dt = D*(d^2*C/d​x^2)-u*(dC​/dx)-R(C). I get the " Index exceeds matrix dimensions" error message, despite properly defining variables.

%% Solving Diffusion-Convection-Reaction PDE (1D) using Crank-Nicholson Method % dC/dt = D*(d^2*C/dx^2)-u*(dC/dx)-R(C) clear all close all clc

%% Parameters xl=0; % Left boundary in space xr=0.14; % Right boundary in space tinitial=0; % Initial time tfinal=30; % Final time tsteps=100; % Number of time steps xsteps=100; % Number of space steps dt=(tfinal-tinitial)/tsteps; % Step sizes in time dx=(xr-xl)/xsteps; % Step sizes in space

%% Coefficients global D u D=3.2e-9; %Diffusion coefficient u=1.6e-2; % Convection velocity

%% Reaction global Stoic k_tot SSA_BET MW_CaCO3 C_CaCO3 C_CaCO3_0 K_ad C_H t Stoic = 1.647; % Stoichiometric coefficient k_tot = 0.01; % Mass transfer coefficient SSA_BET = 12.54; % Specific surface area by BET MW_CaCO3 = 100.0869; % Molecular weight of CaCO3 K_ad = 880.3; % Adsorption constant C_CaCO3_0 = 0.019982635; % Initial concentration of CaCO3 C_CaCO3 = C_CaCO3_0 - (C_CaCO3*C_H*MW_CaCO3*SSA_BET*k_tot*t)/(C_H*K_ad + 1); % Concentration of CaCO3 at t+1 Rxn = lambertw(0, K_ad*exp(log(exp(K_ad/100)) - 2*log(10) - C_CaCO3*MW_CaCO3*SSA_BET*k_tot*Stoic*t))/K_ad; %Rxn = -Stoic*k_tot*SSA_BET*MW_CaCO3*C_CaCO3*C_H*(1-(K_ad*C_H)/(1+K_ad*C_H)); %Rxn = C_H(t,C_H);

%% Concentration of H+ at the boundaries for i=1:tsteps+1 C_H(1,i)=0.; C_H(xsteps+1,i)=0.; t(i)=(i-1)*dt; end

%% Initial concentration of H+ for i=1:xsteps+1 x(i)=(i-1)*dx; C_H(i)=0.01; end

%% Coefficients used to generate matrices global alpha beta gama alpha=dt*D/(2*dx); beta=-dt*u/(4*dx); gama=-dt*Rxn./2;

%% Defining the matrices M-left global aal bbl ccl MMl aar bbr ccr MMr aal(1:xsteps-2)=-alpha-beta; bbl(1:xsteps-1)=1+2.*alpha-gama(1:tsteps-1); ccl(1:xsteps-2)=-alpha+beta; MMl=diag(bbl,0)+diag(aal,-1)+diag(ccl,1);

%% Defining the matrices M-right aar(1:xsteps-2)=alpha+beta; bbr(1:xsteps-1)=1-2.*alpha+gama(1:tsteps-1); ccr(1:xsteps-2)=alpha-beta; MMr=diag(bbl,0)+diag(aal,-1)+diag(ccl,1);

%% Implementation of the Crank-Nicholson method for i=2:tsteps % Time loop C_HC_H=C_H(2:xsteps,i-1); C_H(2:xsteps,i)=inv(MMl)*MMr*C_HC_H; end

%% Reaction %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %function Rxn = C_H(t,C_H) %Stoic = 1.647; %k_tot = 0.01; %SSA_BET = 12.54; %MW_CaCO3 = 100.0869; %K_ad = 880.3; %Rxn = (k_tot*SSA_BET*MW_CaCO3)*(0.019982635-(C_H/Stoic))*C_H*(1-((K_ad*C_H)/(1+K_ad*C_H))); %end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Plotting the results

plot(t,C_H') xlabel('Time (sec)') ylabel('Concentration (mol/dm3)')