Memoryless Weak Nonlinearity Composition

From Wikimization

% Compute composition y(z(x)) coefficients of x given mapping coefficients zeta & xi.
% zeta : memoryless weak power series z(x), first-order leading: [zeta_1, zeta_2 ... zeta_Nt]
% xi   : memoryless weak power series y(z), first-order leading: [  xi_1,   xi_2 ...   xi_Nh]
% Mapping coefficients, constituting nonlinearity composition, zeta & xi are not commutative.
function y = composite2(zeta, xi, precision)
   if nargin < 3, precision = 34; end
   mp.Digits(precision);  %Advanpix Multiprecision Toolbox
   zeta = zeta(:)'; 
   xi   =   xi(:)';
   y       = mp('0');
   z_power = mp('1'); 
   for power = 1:numel(xi)
      z_power = conv(z_power, zeta); 
      term = [zeros(1,power-1,'mp'), xi(power)*z_power];
      y    = term + [y, zeros(1,numel(term)-numel(y),'mp')];
   end
end