function F1= multaylor(f,x,x0,P,w,I)
%MULTAYLOR computes the Taylor expansion of a multivariable vector function
%
% F1 = MULTAYLOR(F,X,X0,P), for vector function F=[F1;...;Fm], where
% Fk are multivariable functions defined in terms of the symbolic variables
% Xk in X=[X1,...,Xm], returns the  m x (m+P)!/(m!P!) coefficient matrix
% of the Pth order Taylor expansion of F around X0.  
%
% F1 = MULTAYLOR(F,X,X0,P,W,I) uses, for faster execution,
% pre-calculated W = MTAYLORWEIGHTS(M,P) and I = MTAYLORIND(M,P)
%
%   See also CARLBILIN  

%   CBVK Toolbox
%   JHMG, PMSB 2016


m = length(x);
if nargin < 5
   w = mtaylorweights(m,P);
end
% Non-redundant indices
if nargin < 6
   I = mtaylorind(m,P);
end

% Numerical variables named after the elements of x and set to x0
for i=1:m
  eval([char(x(i)) '=x0(i);']);
end

MP = nchoosek(m+P,m);
F1 = zeros(m,MP);

S=whos('f');
if strcmp(S.class,'double')
    F1(:,1)=f(:)';
  return
end

fst=char(f);
F1(:,1)=(eval(fst(8:end-1)))';

Js = jacobian(f,x);
for i=1:m

 J = Js(i,:);
  Jst=char(J);
  F1(i,2:m+1)=eval(Jst(8:end-1)); 
 
   ini = m+1;
   for p=2:P
      J = jacobian(J,x);
  
      % Selection of non-redundant terms
      Jr = selnrjac(char(J),I{p-1});
      
      % Evaluation at x0
      Jlin = eval(Jr(8:end-1));
      F1(i,ini+1:ini+length(Jlin))=Jlin;
      ini = ini+length(Jlin);
      
      % Back from string to symbolic
      J=evalin(symengine,Jr);
   end
end

% Weight application
F1 = F1.*repmat(w(1:size(F1,2)),[m 1]);
