Matlab Calling Functions from Mupad Notebooks in MatlabR2015a

After loading "main.m" and the MuPAD notebook, here is the error message I got:

Warning: The system is rank-deficient. Solution is not unique. 
> In symengine (line 57)
  In sym/privBinaryOp (line 903)
  In  /  (line 297)
  In main (line 93) 
Error using lsqfcnchk (line 108)
If FUN is a MATLAB object, it must have an feval method.

Error in lsqnsetup (line 37)
    funfcn = lsqfcnchk(FUN,caller,lengthVarargin,funValCheck,flags.grad);

Error in lsqnonlin (line 182)
[funfcn,mtxmpy,flags,sizes,funValCheck,xstart,lb,ub,EXITFLAG,Resnorm,FVAL,LAMBDA, ...

Error in main (line 101)
X = lsqnonlin(Call(kappa,theta,delta,vega,rho) - Market,X0,[0,0,0,0,0],[],options);

Actually I do not plan to solve functions l(n),m(n),and c(n) explicitly. They are defined recurrently. How should I fix this bug? Can my idea be fulfilled by Matlab? (maybe Matlab would not do the recurrence equations for me?)

Excel is attached to this question.

main.m

% This Task Consists 2 Steps.

% Step1: Debug "main.m" so that for given values of (K,t,s), "Call" can be expressed as a function of 5 parameters:
% (kappa,theta,delta,vega,rho), namely Call = function(kappa,theta,delta,vega,rho). This expression can be implicit.  

% Step2: After reading 100+ sets of numbers (K,t,s) from excel all at once, we get an array of Call's.
% Calibrate the array of Call's to the array of "Market Values" in excel using the Levenberg-Marquardt algorithm,
% a built-in function in Matlab, to get the optimal values of parameters (kappa,theta,delta,vega,rho).

syms kappa theta delta vega rho    % 5 parameters
syms tau n

% Case for 2017.10.4
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

K = xlsread('2017.10.4.xls','B5:B190');
t = xlsread('2017.10.4.xls','C5:C190');
s = 0;
Future = 11.35;
r = 0.0108;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

t0 = 0.00001;           % Fixed data
g = 1/2;                % Fixed data

% Define variables using (kappa,theta,delta,vega,rho) %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
alpha = 2 * kappa / delta^2 + 1;

betaa = 2 * kappa * theta / delta^2 + 1;

C = 1 / sqrt(2 * pi * vega);

lambda = 1 / (2 * vega);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Define function l(n)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms l(n)

% result = feval(symengine,l,n);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Define nested functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% f1 is a function of tau, which disappears in f2
f1 = C * rho * (1 + t0)^(g * rho - 1) * tau^(-g - 1) * exp(-lambda * tau * (1 + t0)^(-rho)) * (g + lambda * tau * (1 + t0)^(-rho));

% f2 depends on f1; f2 is an integral of tau
f2 = int((1 - exp(-n * kappa * theta * tau)) * f1,tau,0,inf);

% f3 depends on f2 and l(n)
f3 = exp(-t * f2) * sqrt(factorial(n) / gamma(n + alpha)) * l(n);

% A depends on f3; n disappears in A
A = betaa * gamma(alpha - 1) * symsum(f3,n,0,inf);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

k = K .* A / Future;

% Define function m(n)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms m(n)

% result = feval(symengine,m,n);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Define function c(n)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms c(n)

% result = feval(symengine,c,n);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Define function B
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B = exp((s - t) * f2) * c(n) * l(n);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Express Call implicitly using (kappa,theta,delta,vega,rho)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X = [kappa,theta,delta,vega,rho];
Call = symfun(exp(-r * t) * Future / A * symsum(B,n,0,inf), X);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Perform the Levenberg-Marquardt algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X0 = [0,0,0,0,0];
Market = xlsread('2017.10.4.xls','F5:F190');
options = optimoptions(@lsqnonlin,'Algorithm','levenberg-marquardt');
X = lsqnonlin(Call(kappa,theta,delta,vega,rho) - Market,X0,[0,0,0,0,0],[],options);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Here is the code in a single MuPAD notebook(hit "Enter" after typing each equation):

solve(rec(l(n) = (alpha - 2 + 2 * n - betaa) / sqrt(n * (alpha - 1 + n)) * l(n - 1) - sqrt((alpha + n - 2) * (n - 1) / (alpha - 1 + n) / n) * l(n - 2),l(n),{l(0) = 1 / sqrt(gamma(alpha)), l(1) = (alpha - betaa) / sqrt(gamma(alpha + 1))}))

solve(rec(m(n) = (alpha + 2 * n - 1 - betaa / k) / sqrt(n * (alpha + n)) * m(n - 1) - sqrt((alpha + n - 1) * (n - 1) / (alpha + n) / n) * m(n - 2),m(n),{m(0) = 1 / sqrt(gamma(alpha + 1)), m(1) = (1 + alpha - betaa / k) / sqrt(gamma(alpha + 2))}))

solve(rec(c(n) = c(n-1) * sqrt(n / (alpha + n - 1)) + betaa^alpha / k^(alpha - 1) * m(n - 2) * exp(-betaa / k) / sqrt(n * (n - 1) * (n + alpha - 1)),c(n),{c(0) = sqrt(gamma(alpha)) * (betaa / (alpha - 1) - k) * gammainc(betaa / k,alpha - 1) + (betaa^(alpha - 1) / k^(alpha - 2)) * exp(-betaa / k) / sqrt(gamma(alpha)),c(1) = sqrt(gamma(alpha + 1)) * betaa / alpha / (alpha - 1) * (1 - igamma(alpha - 1,betaa / k) / gamma(alpha - 1))}))