Using Matlab to calculate numerically the jacobian of a vector-valued function

Question

I would like to use numerical methods to calculate the Jacobian of a vector-valued function in Matlab. Matlab has function for this, called jacobianest

The equation that I am working with requires finding a root. Here is the code

syms cstar; 
vpasolve(1- normcdf(a,mu_t+ ((cstar-mu_c)*sigma_tc)/sigma_c^2,sigma_cond) - (cstar/100)== 0, cstar);
function application_rate  = apply(ParametersIn);
application_rate  = [apply_rate];
apply_rate = normcdf(0,mu_c-cstar,sigma_c);
[jac,err] = jacobianest(@(param) apply(param), [ParametersIn])

I am using the following parameter values:

mu_t       =    ParametersIn(1);
sigma_t    =    ParametersIn(2);
mu_c       =    ParametersIn(3);
sigma_c    =    ParametersIn(4);
sigma_tc   =    ParametersIn(5);
ParametersIn = [45;13;136;192;1700]
sigma_cond  = (0.5*(2*sigma_t^2-((2*sigma_tc^2)/sigma_c^2)))^(1/2);

The jacobianest code is not working – Matlab doesn't seem to be able to deal with the implicit 'cstar' function. This is the error I get:

Error using mupadengine/feval (line 166)
An arithmetical expression is expected.
Error in sym/vpasolve (line 172)
    sol = eng.feval('symobj::vpasolve',eqns,vars,X0);

Any idea how to fix this?