Memoryless Weak Nonlinearity Composition

From Wikimization

% Compute y(z(x)) coefficients of x given mapping coefficients xi and zeta.
% xi   : polynomial y(z), first order leading: [xi_1,   xi_2,   ...]
% zeta : polynomial z(x), first order leading: [zeta_1, zeta_2, ...]
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