function F1 = mtredlit(f,x,P)
%MTREDLIT computes the literal Taylor expansion of a multivariable vector
%function around the origin
%
% F1 = MTREDLIT(F,X,P), for vector function F=[F1;...;Fm], where Fk are
% multivariable literal functions defined in terms of the symbolic
% variables Xk in X=[X1,...,Xm], returns the  m x (1+m+m^2+...+m^P) literal
% coefficient matrix of the Pth order Taylor expansion of F around X0. For
% greater overall efficiency, the factorial constants are left out of the
% coeeficients, only to be included at the numerical evaluation step
%
%   See also MULTAYLORLLIT

%   CBVK Toolbox JHMG, PMSB 2016

warning('cbvk:mtredlit','Consider the use of the multaylorlit, which is possibly faster');
 
m = length(x);

F1 = cell(m,P+1);

fst=evalin(symengine,subsx(char(f)));
for i=1:m
    F1{i,1}=char(fst(i));
end

Js = jacobian(f,x);

for i=1:m
  
  % p=1  
  J = Js(i,:);
  Jst=char(evalin(symengine,subsx(char(J))));
  F1{i,2}=Jst(8:end-1); 

   for p=2:P
     
      J = jacobian(J(:),x);
      
      % Evaluation at the origin
      Jst=char(evalin(symengine,subsx(char(J))));
      F1{i,p+1}=Jst(8:end-1);
   end
end
