C3coeff.m

function C3x = C3coeff(n)
%C3COEFF  Evaluate coefficients for C_3
%
%   C3x = C3COEFF(n) evaluates the coefficients of epsilon^l in Eq. (25).
%   n is a scalar.  C3x is a 1 x 15 array.

  persistent coeff nC3 nC3x
  if isempty(coeff)
    nC3 = 6;
    nC3x = (nC3 * (nC3 - 1)) / 2;
    coeff = [ ...
        3, 128, ...
        2, 5, 128, ...
        -1, 3, 3, 64, ...
        -1, 0, 1, 8, ...
        -1, 1, 4, ...
        5, 256, ...
        1, 3, 128, ...
        -3, -2, 3, 64, ...
        1, -3, 2, 32, ...
        7, 512, ...
        -10, 9, 384, ...
        5, -9, 5, 192, ...
        7, 512, ...
        -14, 7, 512, ...
        21, 2560, ...
            ];
  end
  C3x = zeros(1, nC3x);
  o = 1;
  k = 1;
  for l = 1 : nC3 - 1
    for j = nC3 - 1 : -1 : l
      m = min(nC3 - j - 1, j);
      C3x(k) = polyval(coeff(o : o + m), n) / coeff(o + m + 1);
      k = k + 1;
      o = o + m + 2;
    end
  end
end