function [F1i,F1c,F11] = mtredsp(f,x,P)
%MTREDSP computes the literal Taylor expansion of a multivariable vector
%function around the origin. The expansion is returned in a sparsity
%oriented format
%
% [F1I,F1C,F11] = MTREDSP(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. F11
% is a string with the first column of coeffcient matrix. F1i and F1C are a
% sparsity oriented cell representation of the remaining columns
%
%   See also MULTAYLORSP  

%   CBVK Toolbox
%   JHMG, PMSB 2016

warning('cbvk:mtredsp','Consider the use of the multaylorsp, which is possibly faster');

m = length(f);
mp=cumsum(m.^(1:P));

F1i = cell(1,m);
F1c=F1i;

fst=char(evalin(symengine,subsx(char(f))));
F11=fst(8:end-1);

% Jacobian in symbolic format
Jsy = jacobian(f,x);

for i=1:m
  J = Jsy(i,:);
  
  % p=1
  % Evaluation at the origin
  Jst=char(evalin(symengine,subsx(char(J))));
  F1c{i}='[';
  F1i{i}=[];

  
  % Selection of non-null terms
  [fct,fit]=spar(Jst(8:end-1),[0 mp(1)]);
  if ~isempty(fit)
    F1i{i}=[F1i{i} fit];
    F1c{i}=[F1c{i}  fct];
  end
  
  
  for p=2:P
     
      J = jacobian(J(:),x);
    
    Jst=char(evalin(symengine,subsx(char(J))));
    
    % Sparsity
    [fct,fit]=spar(Jst(8:end-1),[mp(p-1) mp(p)]);
    if ~isempty(fit)
      F1i{i}=[F1i{i} fit];
      F1c{i}=[F1c{i} ',' sprintf('%0.15g',1/factorial(p)) '*' fct];
    end
  end
   F1c{i}=[F1c{i} ']'];
end