function F1= multaylorlit(f,x,P,I)
%MULTAYLORLIT computes the literal Taylor expansion of a multivariable
%vector function around the origin
%
% F1 = MULTAYLORLIT(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 in cell form the  blocks of the m x
% (m+P)!/(m!P!) literal coefficient matrix of the Pth order Taylor
% expansion of F around the origin.
%
% F1 = MULTAYLORLIT(F,X,P,I) uses, for faster execution, pre-calculated I =
% MTAYLORIND(M,P)
%
%   See also CARLBILIN
%
%   CBVK Toolbox JHMG, PMSB 2016


m = length(x);
if nargin < 4
   I = mtaylorind(m,P);
end

F1 = cell(m,P+1);
fst=evalin(symengine,subsx(char(f)));

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

Jsy = jacobian(f,x);
for i=1:m
  
  % p=1  
  J = Jsy(i,:);
  Jst=char(evalin(symengine,subsx(char(J))));
  F1{i,2}=Jst(8:end-1);
  
  for p=2:P
     
    J = jacobian(J(:),x);
    
    % Selection of non-redundant terms
    Jr = selnrjac(char(J),I{p-1});
    
    % Back from string to symbolic
    J=evalin(symengine,Jr);
    
    % Evaluation at the origin
    Jst=char(evalin(symengine,subsx(Jr)));
    F1{i,p+1}=Jst(8:end-1);
    
  end
end
