ODE23s Passing Additional Arguments to Symbolic Jacobian

I am passing a symbolic Jacobian to an ODE23s anon fun and I am not entirely sure how to include additional vars which need to be substituted into the Jacobian in the intialization of the Jacobian. I do not have the entire code or function as it is a subsection of working code.

The call back is feval step, yet there is no option to pass additional arguments to the Jacobian.

QUESTION: How do I pass additional args to symbolic vars in the Jacobian of an ODE solver?

(I have attempted to edit and run my own ODE file but this seems like it skips a step or 2.)

Error message:

Not enough input arguments.

Error in symengine>@(t,y,L_f,theta_h,TsubsODE,ki,kw,q2_ode,qloss_ode,htm1,mdot_net_ode,Btm1,rho_w,rho_i,delta_t)[-(kw.*1.0./y.^2.*(TsubsODE-theta_h)-(ki.*qloss_ode.*rho_w.*(q2_ode-qloss_ode.*theta_h).*1.0./(ki+qloss_ode.*(Btm1+(delta_t.*mdot_net_ode)./rho_i+(rho_w.*(htm1-y))./rho_i)).^2)./rho_i)./(L_f.*rho_w),0.0]

Error in ode23s (line 220)
  dfdy = feval(Jac,t,y,Jargs{:});

Error in EMM_C_IceCrystalsRunningWet (line 384)
   [~,b] = ode23s(@(t,b) ode45_EMM_fcn_dhdt(t,y,L_f,theta_h(pRW),TsubsODE,ki,kw,q2_ode,qloss_ode,htm1,mdot_net_ode,Btm1,rho_w,rho_i,delta_t) , [t t+delta_t] , p0_all, opts2); % internal water layer thickness. ode23s( anonmyous function of t,y  function , time span for integration  ,  initial condition (y0) at t0)

Code:

opts = odeset('RelTol',1e-4, 'AbsTol',5e-4, 'Jacobian', JacobSymEMM(L_f,theta_h(pRW),TsubsODE,ki,kw,q2_ode,qloss_ode,htm1,mdot_net_ode,Btm1,rho_w,rho_i,delta_t), 'Vectorized', 'on');

[~,y] = ode23s(@(t,y) ode45_EMM_fcn_dhdt(t,y,L_f,theta_h,TsubsODE,ki,kw,q2_ode,qloss_ode,htm1,mdot_net_ode,Btm1,rho_w,rho_i,delta_t) , [t t+delta_t] , p0_all, opts); % internal water layer thickness. ode23s( anonmyous function of t,y  function , time span for integration  ,  initial condition (y0) at t0)

function Jacob = JacobSymEMM(t,y,L_f,theta_h,TsubsODE,ki,kw,q2_ode,qloss_ode,htm1,mdot_net_ode,Btm1,rho_w,rho_i,delta_t)
syms t y L_f theta_h TsubsODE ki kw q2_ode qloss_ode htm1 mdot_net_ode Btm1 rho_w rho_i delta_t
dhdt = (1/(rho_w.*L_f)).* (  -(kw./y).*(theta_h-TsubsODE) +  ki.*(q2_ode - theta_h.*qloss_ode)./(ki + qloss_ode.*(Btm1 + (1/rho_i).*mdot_net_ode.*delta_t + (rho_w./rho_i).*(htm1 - y))));
JacobianAnal = jacobian(dhdt,[y, t]);
Jacob = matlabFunction(JacobianAnal, 'vars', [t,y, L_f,theta_h,TsubsODE,ki,kw,q2_ode,qloss_ode,htm1,mdot_net_ode,Btm1,rho_w,rho_i,delta_t]);
end



L_f =

334000


theta_h =

273.1500


TsubsODE =

283.1500


ki =

2.2200


kw =

0.6080


q2_ode =

2.0887e+05


qloss_ode =

1.0568e+03


htm1 =

0


mdot_net_ode =

0.2525


Btm1 =

0


rho_w =

1000


rho_i =

920


delta_t =

0.5000