I am trying to solve a combination of algebraic and differential equations . But I am getting an error of "Variables must be specified by univariate function calls"

My MATLAB version is 2018a. I mimicked the code from the MATLAB topic "Solve Differential Algebraic Equations (DAEs)". But i don't know why I am getting this error of "Variables must be specified by univariate function calls".

syms n0(t) n1(t) n2(t) n3(t) p0(t) p1(t) p2(t) p3(t) s0(t) s1(t) s2(t) vac2n(t) vac2p(t) vac2s(t) c2atx0(t) c2atxL(t) dn1(t) dn2(t) dp1(t) dp2(t) vac2n vac2s vac2p c2atx0 c2atxL

syms eta2n eta2p eta2s D2n D2p D2s eta1n eta1p rn rp an ap tplus F c20 I T R vajn vajp k2s k2p k2n

syms csnmax cspmax csn0 csp0 kn kp jn0 jp0 sigma1n sigma1p ln lp lls jnatx0 jpatxL t

eqn1 = n1(t) == 0;

eqn2 = p1(t) == 0;

eqn3 = s0(t) == n0(t);

eqn4 = p0(t) == s2(t)*(lls^2) + s1(t)*lls + s0(t);

eqn5 = D2n*(-2*n2(t)*ln-3*n3(t)*(ln^2)) == D2s*(2*s2(t)*ln+s1(t));

eqn6 = D2n*(2*p2(t)*lp+3*p3(t)*(lp^2)) == D2s*(2*s2(t)*(ln+lls)+s1(t));

eqn7 = vac2n == n0(t) + (2*n2(t)*(ln^2))/3 + (3*n3(t)*(ln^3))/4;

eqn8 = vac2p == p0(t) + (2*p2(t)*(lp^2))/3 + (3*p3(t)*(lp^3))/4;

eqn9 = vac2s == s0(t) + (2*s2(t)*(lls^2))/3 + (s1(t)*lls)/2;

eqn10 = ln*eta2n*vac2n + lp*eta2p*vac2p + lls*eta2s*vac2s == (ln*eta2n+lp*eta2p+lls*eta2s)*c20;

eqn11 = c2atx0 == n2(t)*(ln^2) + n3(t)*(ln^3) + n0(t);

eqn12 = c2atxL == p2(t)*(lp^2) + p3(t)*(lp^3) + p0(t);

eqn13 = eta2n*diff(vac2n,t) == -D2n*(2*n2(t)+3*n3(t)*ln)+an*(1-tplus)*vajn;

eqn14 = eta2p*diff(vac2p,t) == -D2p*(2*p2(t)+3*p3(t)*lp)+ap*(1-tplus)*vajp;

eqn15 = eta2n*diff(c2atx0,t) == -2*D2n*n2(t)+an*(1-tplus)*jnatx0;

eqn16 = eta2p*diff(c2atxL,t) == -2*D2p*p2(t)+ap*(1-tplus)*jpatxL;

eqns = [eqn1 eqn2 eqn3 eqn4 eqn5 eqn6 eqn7 eqn8 eqn9 eqn10 eqn11 eqn12 eqn13 eqn14 eqn15 eqn16];

vars = [n0(t); n1(t); n2(t); n3(t); p0(t); p1(t); p2(t); p3(t); s0(t); s1(t); s2(t); vac2n; vac2p; vac2s; c2atx0; c2atxL];

origVars = length(vars);

M = incidenceMatrix(eqns,vars)

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

% isLowIndexDAE(eqns,vars)

% [DAEs,DAEvars] = reduceDAEIndex(eqns,vars)

% [DAEs,DAEvars] = reduceRedundancies(DAEs,DAEvars)

peqns = symvar(eqns);

peqnvars = symvar(vars);

extraParams = setdiff(peqns,peqnvars)

f = daeFunction(eqns,vars,extraParams); %% I am getting error in this line.

The error is as follows:

Error using mupadengine/feval (line 187)

Variables must be specified by univariate function calls.

Error in sym/daeFunction (line 117)

A = feval(symengine, 'daetools::daeFunction', eqs, vars, params{:});

Error in algebraicsolve (line 39)

f = daeFunction(eqns,vars,extraParams);