function [fi,fc,F11] = multaylorsp(f,x,P,I)
%MULTAYLORSP computes the literal Taylor expansion of a multivariable
%vector function around the origin. The expansion is returned in a sparsity
%oriented format.
%
% [fi,fc,F11]  = MULTAYLORSP(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  blocks of the m x
% (m+P)!/(m!P!) literal coefficient matrix of the Pth order Taylor
% expansion of F around the origin. F11 is the first column of coeffcient
% matrix. Fi and FC are a sparsity oriented representation of the remaining
% columns.
%
% [fi,fc,F11]  = MULTAYLORSP(F,X,P,I) uses, for faster execution,
% pre-calculated I = MTAYLORIND(M,P)

%   CBVK Toolbox JHMG, PMSB 2016

m = length(f);
mp=cumsum(fliplr(calcnumterms(P,m)));

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


fi = cell(1,m);
fc=fi;

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))));
  fc{i}='[';
  fi{i}=[];
  
  % Selection of non-null terms
  [fct,fit]=spar(Jst(8:end-1),[0 mp(1)]);
  if ~isempty(fit)
    fi{i}=[fi{i} fit];
    fc{i}=[fc{i}  fct];
  end
  
  
  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)));
  
    % Sparsity
    [fct,fit]=spar(Jst(8:end-1),[mp(p-1) mp(p)]);
    if ~isempty(fit)
      fi{i}=[fi{i} fit];
      fc{i}=[fc{i} ',' fct];
    end
  end
   fc{i}=[fc{i} ']'];
end
