Index in position 1 exceeds array bounds (must not exceed 1)

Hi everyone, when I run the code below (or in the attachment), I get the error message:

"

Index in position 1 exceeds array bounds (must not exceed 1).

Error in sym/subsref (line 859)
            R_tilde = builtin('subsref',L_tilde,Idx);

Error in
symengine>@(x,in2,in3,param1,param2,param3,param4,param5,param6,param7,param8,param9,param10,param11,param12,param13)[in3(6,:).*param1+in3(7,:).*param2+in3(8,:).*param3;-in3(7,:).*param2-in3(8,:).*param3.*2.0+in3(9,:).*param4-in3(10,:).*param5;param6-(in2(2,:).*in2(4,:))./in2(1,:);param7-(in2(3,:).*in2(4,:))./in2(2,:);param8-in2(4,:).*in2(5,:);in2(6,:)-in3(1,:);in2(7,:)-in3(2,:);in2(8,:)-in3(3,:);in2(9,:)-in3(4,:);in2(10,:)-in3(5,:)]

Error in samplecode>@(t,Y,YP)f(t,Y,YP,A,B,C,D,E,FF,G,H,I,J,K,LL,MM)

Error in sym>funchandle2ref (line 1291)
    S = x(S{:});

Error in sym>tomupad (line 1204)
    x = funchandle2ref(x);

Error in sym (line 211)
                S.s = tomupad(x);

Error in checkDAEInput (line 25)
   eqs = sym(eqs);

Error in sym/decic (line 144)
[eqs, vars] = checkDAEInput(eqs, vars);

Error in samplecode (line 59)
[y0, yp0] = decic(F, 0, y0est, [], yp0est, [], opt);
"

when I debbug, it point me here:

"

[y0, yp0] = decic(F, 0, y0est, [], yp0est, [], opt);

"

However I don't see anything wrong with this entry. What can I do to fix this error?

The complete code is:

clear all
close all
clc

syms a(x) b(x) c(x) d(x) e(x) A B C D E FF G H I J K LL MM

eqn1 = A * diff(a(x),x,2) + B * diff(b(x),x,2) + C * diff(c(x),x,2) == 0;

eqn2 = D * diff(d(x),x,2) - B * diff(b(x),x,2) - 2 * C * diff(c(x),x,2) - E * diff(e(x),x,2) == 0;

eqn3 = FF == (d(x) * b(x)) / a(x);

eqn4 = G == (d(x) * c(x)) / b(x);

eqn5 = H == d(x) * e(x);

eqns = [eqn1 eqn2 eqn3 eqn4 eqn5];

vars = [a(x); b(x); c(x); d(x); e(x)];

origVars = length(vars);

M = incidenceMatrix(eqns, vars);

[eqns, vars] = reduceDifferentialOrder(eqns, vars);

isLowIndexDAE(eqns,vars);

f = daeFunction(eqns,vars, A, B, C, D, E, FF, G, H, I, J, K, LL, MM);

A = 2.89e-9;   
B = 2.35e-9;  
C = 1.69e-9;  
D = 1.6e-8;    
E = 9.25e-9; 
FF = 6.24;      
G = 5.68e-5;  
H = 5.3e-8; 
I = 4.01e-5; 
J = 7.21e-05 ; 
K = 149;
LL = 84.9;
MM = 3.47;

F = @(t, Y, YP) f(t, Y, YP, A, B, C, D, E, FF, G, H, I, J, K, LL, MM);

vars;

a0 = K * LL;
b0 = (FF * K * LL)/MM;
c0 = (FF * G * K * LL)/(MM)^2;
d0 = sym(3.47);
e0 = H/3.47;

y0est = [a0; b0; c0; d0; e0];
yp0est = zeros(5,1);
opt = odeset('RelTol', 10.0^(-2), 'AbsTol' , 10.0^(-2));
[y0, yp0] = decic(F, 0, y0est, [], yp0est, [], opt);
[tSol,ySol] = ode15i(F, [0, 7.21e-05], y0, yp0, opt);
for k = 1:origVars
  S{k} = char(vars(k));
end

plot(xSol,ySol(:,1),'r:',xSol,ySol(:,2),'k-.',xSol,ySol(:,3),'b--', xSol,ySol(:,4),'c-',xSol,ySol(:,5),'m-')
ylabel('Conc (mol/m^{3}')
xlabel ('Time (sec)')
legend('SO_2','HSO_{3}^{-}', 'SO_{3}^{2-}', 'H^+','OH^-','location','Best')
clear all