Problem with "index exceeds array bound"

Hello! I'm trying to model turbulent mixing with chemical reactions using engulfment model and ode45. To do this I need to discretise into sigma parts the added volume of acid to the solution. I must solve a system of ordinary differential equations for each part of the feed (each time until one of reagents is consumed). From one part to another, I need to update concentrations of reagents in the vicinity of the reaction zone (for the (k+1)th part of the feed concentrations are Ci(k), but I changed matrix indices, because I can't start from 0). I'd like to save all the final values of variables for each part. These values are essential to calculate Ci(k+1) for the next part of the feed.

The script usually stops after calculation for the 1st part. The error which constantly occurs is "index exceeds array bounds" (additionally with errors concerning used functions). I've tried to change indexing, but without success. I currently have no idea how to properly fix this, so I'd greatly appreciate any remarks and suggestions. Thank you!

This is my event function:

function [value,isterminal,direction] = eventsfunction(~,c,~)
value = c(3);
isterminal = 1;
direction = -1;
end

This is main function:

function [dcdt] = mojafunkcja(t,c,k)
% standard gravity [m/s^2]
g = 9.80665;
% impeller diameter [m]
Dm = 0.052;
% average temperature [K]
T = 273.15+20; 
% density [kg/m^3]
ro = 998.2063; 
% dynamic visosity [Pa*s]
M1 = -52.843;
M2 = 3703.6;
M3 = 5.866;
M4 = -5.879*10^(-29);
M5 = 10;
mi = exp(M1+M2/T+M3*log(T)+M4*T^M5);
% kinematic viscosity [m^2s]
ni = mi/ro;
% impeller speed [1/s]
Nm = 500/60;
% Reynolds number
Rem = (Nm*Dm^2)/ni;
% Froude number
Frm = (Dm*Nm^2)/g;
% Euler number
Eum = 1.1*Frm^((1-log10(Rem))/40);
% initial volume of solution [dm^3]
Vo = 0.2;
% average energy dissipation rate [W/kg]
epsilon = (Eum*Nm^3*Dm^5)/(Vo*10^(-3));
% added volume of acid solution [dm^3]
Vbo = 0.001;
% E parameter
E = 0.05776*(epsilon/ni)^(1/2);
% inital ionic force of solution [mol/dm^3]
Io = 0.10740;
% rate constant for Dushman reaction [(dm^3)^4/(mol^4*s)]
k2 = 10^(2.1766*Io+8.9484);
% rate constant for 3rd reaction from left to right [dm^3/(mol*s)]
k3 = 5.6*10^9;
% rate constant for 3rd reaction from right to left [1/s)]
k33 = 7.5*10^6;
% feed time [s]
tf = 4.5;
% number of parts
sigma = 100;
% volumetric flow [m^3/s]
Q = (Vbo*0.001/sigma)/tf;
% average fluid velocity in the bulk [m/s]
u = pi*Dm*Nm;
% integral scale for fluctuations of concentrations [m]
Lambdac = (Q/(pi*u))^(1/2);
% mesomixing time [s]
ts = 2*(Lambdac^2/epsilon)^(1/3);
% initial molar concentrations of reagents in reaction zone [mol/dm^3]  
C1(1) = 0.01179;
C2(1) = 0.00234;
C3(1) = 0.0933;
C4(1) = 0.0898;
C5(1) = 0;
C6(1) = 0;
dcdt = zeros(8,1);
% 1-I-,2-IO3-,3-H+(C3-H2BO3-),4-H3BO3,5-I2,6-I3-,7-Vb,8-Xb
dcdt(1) = E*(C1(k)-c(1))*(1-(c(7)*exp(-t/ts))/(Vbo/sigma))-5*k2*(c(1)^2)*c(2)*(c(3)^2)-k3*c(5)*c(1)+k33*c(6);   
dcdt(2) = E*(C2(k)-c(2))*(1-(c(7)*exp(-t/ts))/(Vbo/sigma))-k2*(c(1)^2)*c(2)*(c(3)^2);
dcdt(3) = -E*C3(k)*(1-(c(7)*exp(-t/ts))/(Vbo/sigma))-6*k2*(c(1)^2)*c(2)*(c(3)^2);
dcdt(4) = E*C3(k)*(1-(c(7)*exp(-t/ts))/(Vbo/sigma))+E*(C4(k)-c(4))*(1-(c(7)*exp(-t/ts)));
dcdt(5) = E*(C5(k)-c(5))*(1-(c(7)*exp(-t/ts))/(Vbo/sigma))+3*k2*(c(1)^2)*c(2)*(c(3)^2)-k3*c(5)*c(1)+k33*c(6);
dcdt(6) = E*(C6(k)-c(6))*(1-(c(7)*exp(-t/ts))/(Vbo/sigma))+k3*c(5)*c(1)-k33*c(6);
dcdt(7) = E*c(7)*(1-(c(7)*exp(-t/ts))/(Vbo/sigma));
dcdt(8) = E*c(8)*(1-(c(7)*exp(-t/ts))/(Vbo/sigma));
end

This is my script:

clc;
% initial volume of solution [dm^3]
Vo = 0.2;
% added volume of acid solution [dm^3]
Vbo = 0.001;
% number of parts
sigma = 100;
% time step [s]
h = 10^(-6);
% initial time [s]
t0 = 0;
% time span [s]
deltat = 1;
% final time [s]
tk = t0+deltat;
% 1-I-,2-IO3-,3-H+(C3-H2BO3-),4-H3BO3,5-I2,6-I3-,7-Vb,8-Xb
% initial molar concentrations of reagents in vicinity of reaction zone 
% [mol/dm^3]
c10 = 0;
c20 = 0;
c30 = 0.41273;
c40 = 0;
c50 = 0;
c60 = 0;
% initial volume of reaction zone
c70 = Vbo/sigma;
% initial volume fraction of reaction zone
c80 = (Vbo/sigma)/(Vo+Vbo/sigma);
% initial conditions
c0 = [c10 c20 c30 c40 c50 c60 c70 c80];
% time span
tspan = t0:h:tk;
for k = 1:(sigma)
% initial molar concentrations of reagents in reaction zone [mol/dm^3]  
C1(1) = 0.01179;
C2(1) = 0.00234;
C3(1) = 0.0933;
C4(1) = 0.0898;
C5(1) = 0;
C6(1) = 0;
% c(3)=0
options = odeset('Events',@eventsfunction);
% used function
[t,c,tk,ck] = ode45(@(t,c) mojafunkcja(t,c,k),tspan,c0,options);
tkk(k) = tk;
% final molar concentrations in reaction zone [mol/dm^3]
c1k(k) = ck(1,1);
c2k(k) = ck(1,2);
c3k(k) = ck(1,3);
c4k(k) = ck(1,4);
c5k(k) = ck(1,5);
c6k(k) = ck(1,6);
% final volume of reaction zone [dm^3]
c7k(k) = ck(1,7);
% final volume fraction of reaction zone
c8k(k) = ck(1,8);
% updated molar concentrations of reagents in vicinity of reaction zone
% [mol/dm^3]
C1(k+1) = C1(k)*(1-c8k(k))+c1k(k)*c8k(k);
C2(k+1) = C2(k)*(1-c8k(k))+c2k(k)*c8k(k);
C3(k+1) = C3(k)*(1-c8k(k));
C4(k+1) = C4(k)*(1-c8k(k))+c4k(k)*c8k(k);
C5(k+1) = C5(k)*(1-c8k(k))+c5k(k)*c8k(k);
C6(k+1) = C6(k)*(1-c8k(k))+c6k(k)*c8k(k);
% total volume of solution in vessel [dm^3]
Vc(k) = Vo+(Vbo*k)/(sigma);
% segregation index
Xs(k) = ((C6(k+1)+C5(k+1))*Vc(k))/(c30*(Vbo*k/(sigma)))*(2+C3(1)/(3*C2(1)));        
end

[Adam]

Learn to debug your code, in particular using the 'Pause on Errors' (or 'Stop on Errors' in older versions of Matlab) option from the Breakpoints menu. It will stop on the line that causes the error then you can trivially check on the command line and in the workspace what the sizes of your components and the indices into them are at that point.

It will save you so much time with these type of bugs in future if you are able to debug.

and for now, if nothing else, it will give you something to try until someone else may spot the exact problem.