Error of function (eps: class must be double) trying to solve a symbolic DAE system with ODE15i

Many greetings to anyone who stumbles upon this thread. I'm currently trying to solve thisalgebraic-differential equations system for a chemical reaction rate regression analysis. [OH] is a variable x(t) and [DCF] is my variable y(t). All the other terms are simple constants.

And to do so I stumbled upon this guide for solving DAEs symbolically: https://www.mathworks.com/help/symbolic/solve-differential-algebraic-equations.html

For those that have tried this approach I'd like to ask for help. During Step 4 where I'm instructed to find optimal initial conditions for my variables I get the following errors:

Error using eps
Class must be 'single' or 'double'.

Error in odenumjac (line 88)
br = eps(classF) ^ (0.875);

Error in ode15ipdupdate (line 25)
  [dfdy,dfdy_options.fac,NF] = odenumjac(odefun,{t,y,yp,extras{:}},f,dfdy_options);

Error in ode15ipdinit (line 76)
      ode15ipdupdate([],odefun,t0,y0,yp0,f0,dfdy,dfdyp,dfdy_options,dfdyp_options,extras);

Error in decic (line 70)
    ode15ipdinit(odefun,t0,y0,yp0,res,options,varargin);

Error in parametrows (line 63)
[y0,yp0] = decic(F,0,y0est,[],yp0est,[])

My code is pretty much a copy and paste of that guide, except with my own equations. I honestly don't know how to debug or fix it. I get the same error of the eps function even if I don't try to find optimal conditions and rather jump straight to ode15i. Any idea why? Here's my code:

datos %loading of data
error_residual %loading of other data
%here I begin
syms x(t) y(t) FeOOH OZ ka kb kc kd ke kf kg kh ki kj kk kl k_o k_oh k_nine
eqn1 = ka*x(t)^3+x(t)^2*(kb*y(t)+kc*FeOOH+kd*FeOOH/y(t)-ke*OZ)+x(t)*(kf*OZ*y(t)+kg*(OZ^2)-kh*OZ*FeOOH+ki*(OZ*FeOOH/y(t)))+(-kj*OZ^2*y(t)-kk*OZ^2*FeOOH+kl*(OZ^2*FeOOH/y(t))) == 0;
eqn2 = diff(y(t),1) == -k_o*OZ*y(t)-k_oh*y(t)*x(t)-k_nine*FeOOH*OZ;
eqns = [eqn1 eqn2];
vars = [x(t);y(t)];
origVars = length(vars);
M = incidenceMatrix(eqns,vars);
[DAEs, DAEvars] = reduceDAEIndex(eqns,vars);
[DAEs,DAEvars] = reduceRedundancies(DAEs,DAEvars);
pDAEs = symvar(DAEs);
pDAEvars = symvar(DAEvars);
extraParams = setdiff(pDAEs,pDAEvars);
f = daeFunction(DAEs,DAEvars,FeOOH,OZ,k_o,k_nine,k_oh,ka,kb,kc,kd,ke,kf,kg,kh,ki,kj,kk,kl);
k_o = double(x_parameters(1)); %x_parameters is the array with the data for the constants
k_oh = double(x_parameters(2));
k_9 = double(x_parameters(3));
ka =double(x_parameters(4));
kb = double(x_parameters(5));
kc =double(x_parameters(6));
kd =double(x_parameters(7));
ke = double(x_parameters(8));
kf =double(x_parameters(9));
kg = double(x_parameters(10));
kh = double(x_parameters(11));
ki = double(x_parameters(12));
kj = double(x_parameters(13));
kk = double(x_parameters(14));
kl = double(x_parameters(15));
FeOOH = double(g(1)); %likewise with x_parameters, g is data array
OZ = double(ox(1));

F = @(t,Y,YP) f(t,Y,YP,FeOOH,OZ,k_o,k_nine,k_oh,ka,kb,kc,kd,ke,kf,kg,kh,ki,kj,kk,kl);
y0est = [0;0.054768223996110]; %this numbers are my initial numerical guesses for the initial conditions
yp0est = [0;-0.012990851139745];
[y0,yp0] = decic(F,0,y0est,[],yp0est,[])

I literally followed that guide step by step, as well as the one for solving DAEs using Mass matrix: https://www.mathworks.com/help/symbolic/solve-daes-using-mass-matrix-solvers.html. I get the same error with this one.

[Walter Roberson]

Please attach the two files so that we can run your code

[Miguel Figueroa]

It's a bit messy but here it is, don't pay mind to datos.m and error_residual.m. It's not ver well optimized. Just run it and witness the error of the main script.

edit:s sent the wrong one, fixed.

[Miguel Figueroa]

Hello sir, I've seen to notice what my mistake was, I had specified one symbolic constant as k_nine and later on gave it a numeric value as k_9 = constant. So with that fixed I get the following even more crushing error:

Edit: I'm attaching the modified codes that give this error.

Error using decic (line 108)
Convergence failure in DECIC.

Error in parametrows (line 76)
[y0,yp0] = decic(F,0,double(y0est),[],double(yp0est),[])