Index exceeds array bounds error message in a three-point bvp

Hi everyone, I'm trying to solve a three-point bvp, following:

https://www.mathworks.com/help/releases/R2018a/matlab/ref/bvpinit.html

However, I get the error message:

"Index exceeds array bounds.

Error in ConcentrationProfilesinLimestoneSlurryLiquidFilm>materialbalanceodes (line 63) dydx(1) = y(6);

Error in bvparguments (line 105) testODE = ode(x1,y1,odeExtras{:});

Error in bvp4c (line 130) bvparguments(solver_name,ode,bc,solinit,options,varargin);

Error in ConcentrationProfilesinLimestoneSlurryLiquidFilm (line 4) sol = bvp4c(@materialbalanceodes,@materialbalancebc,solinit);"

What can I try? Please see the code on the attachment or below:

x = [0, 10e-6, 20e-6, 20e-6, 30e-6, 40e-6];

solinit = bvpinit(x,[4.03e-08, 3.39e-02, 2.60e-01, 5.45e-01, 4.82e-03, 0, 0, 0]); sol = bvp4c(@materialbalanceodes,@materialbalancebc,solinit);

y = deval(sol,x); plot(x,y(1,:),x,y(2,:),x,y(3,:),x,y(4,:),x,y(5,:)) legend('SO_2','HSO_{3}^{-}','SO_{3}^{2-}','HCO_{3}^{-}','CO_{3}^{2-}')

global A B C D E F G H I J K L M N O P Q R S T U

%PARAMETERS A = sym(9.4e-4); % ks = solid-liquid mass transfer (m/s) B = sym(3.69e2); % Ap = solid-liquid surface area (m2/m3) C = sym(2.89e-9); % DSO2 = diffusivity of SO2 (m2/s) D = sym(1.6e-9); % DCO32 = diffusivity of CO32- (m2/s) E = sym(4.92e-5); % CBs = concentration of Ca2+ on limestone surface (mol/m3) F = sym(1.6e-8); % DH = diffusivity of H+ (m2/s) G = sym(6.24); % KSO2 = SO2 dissociation equilibrium (mol/m3) H = sym(2.35e-9); % DHSO3 = diffusivity of HSO3 (m2/s) I = sym(5.68e-5); % KHSO3 = HSO3 dissociation equilibrium constant (mol/m3) J = sym(1.69e-9); % DSO32 = diffusivity of SO32 (m2/s) K = sym(149); % HSO2 = Henry's constant (m3Pa/mol) L = sym(47.28325); % PSO2 = partial pressure of SO2 (Pa) M = sym(4.14e-6); % kga = gas-side mass transfer product (1/s) N = sym(8.314); % R = universal gas constant (J/mol.K) O = sym(323.15); % T = reactor tempeature (K) P = sym(2.09); % CH = concentration of H+ at the gas-liquid interface (mol/m3) Q = sym(9.333E-06); % CH = concentration of H+ at the in the bulk slurry (mol/m3) R = sym(0.679); % S(IV) = total sulfur concentration in the bulk slurry (mol/m3) S = sym(4.88e-6); % C(IV) = total sulfur concentration in the bulk slurry (mol/m3) T = sym(1.7e-3); % KCO2 = CO2 dissociation equilibrium constant (mol/m3) U = sym(6.55e-8); % KHCO3 = HCO3 dissociation equilibrium constant (mol/m3)

% VARIABLES

% y(1) = CSO2 % y(2) = CHSO3 % y(3) = CSO32 % y(4) = CHCO32 % y(5) = CCO32

% EQUATIONS

% Region I

% dy(1)/dx = y(6); % d^2y(1)/dx^2 = (2*A*B / C)*(1 + (C*y(1) / 2*D*E) + (F*G*(y(2)/y(1)) / D*E)*E; % dy(2)/dx = y(7); % d^2y(2)/dx^2 = -(2*A*B / C)*(1 + (C*y(1) / 2*D*E) + (F*G*(y(2)/y(1)) /D*E)*E;

%Region II

% dy(3)/dx = y(8) % dy(4)/dx = y(4); % d^2y(3)/dx^2 = -(A*B / D)*(1 + (F*I*(y(3) / y(2)))/(D*E))*E; % dy(5)/dx = y(5)

function dydx = materialbalanceodes(y,x, A ,B, C, D, E, F, G, H, I)

     dydx = zeros(8,1);
      dydx(1) = y(6);
      dydx(2) = (2 * A * B / C) * (1 + ((C * y(1)) / (2 * D * E)) + ((F * G )* (y(2)/y(1)) / (D*E)))*E;
      dydx(3) = y(7);
      dydx(4) = -((2 * A * B) / (H)) * (1 + ((C * y(1)) / (2 * D * E)) + ((F * G) * (y(2) / y(1)) / (D * E))) * E;
      dydx(5) = y(8);
      dydx(6) = y(4);
      dydx(7) = -((A * B) / D) * (1 + (F * I * (y(3) / y(2)))/(D * E)) * E;
      dydx(8) = y(5);  
     end

function res = materialbalancebc(ya,yb, J, K, L, M, N, O, P, Q) res = [(ya(1)- (K * L)); (ya(2)-((I * K * L) / P)); (ya(6)+ M * ((L /( N * O))- (K * L))); (yb(2) - ((R * G * Q)/(Q.^2 + G * Q + G * I))); (yb(3) - ((R * G * I)/(Q.^2 + G * Q + G * I))) ; yb(4) - ((S * T * Q)/(Q.^2 + T * Q + T * U))]; end