function F1= mtred(f,x,x0,P)
%MTRED computes the Taylor expansion of a multivariable vector function
%
% F1 = MTRED(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 (1+m+m^2+...+m^P) coefficient matrix of
% the Pth order Taylor expansion of F around X0.
%
%   See also MULTAYLOR

%   CBVK Toolbox JHMG, PMSB 2016

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


m = length(x);

% for i=1:m
%   varname = sprintf('x%d',i);
%   eval([varname '=x0(i);']);
% end

for i=1:m
  eval([char(x(i)) '=x0(i);']);
end

MP = sum(m.^(0:P));
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);

mps = cumsum(m.^(0:P))+1;


for i=1:m
  J = Js(i,:);
  Jst=char(J);
  F1(i,2:m+1)=eval(Jst(8:end-1)); 

   for p=2:P
     
      p1=1/factorial(p);
      J = jacobian(J(:),x);
      Jst=char(J);
      F1(i,mps(p):mps(p+1)-1)=p1*eval(Jst(8:end-1));
      
   end
end
