Calculating the product of discrete transfer functions I'm trying to construct a higher order IIR from a biquad cascade, I know this is generally frowned upon, but it's necessary as the coefficients will be embedded into a larger formula used in a FDTD boundary condition.

According to this answer, I simply need to multiply the biquad transfer functions together, however I'm not sure how to do that.  Is it simply a case of multiplying each numerator polynomial into a higher order one, and then doing the same for the denominator?
The total transfer function is indeed just the product of the individual transfer functions of the biquads:

$$H(z)=\prod_{i=1}^NH_i(z)\tag{1}$$

If $H(z)=B(z)/A(z)$, this means that

$$B(z)=\prod_{i=1}^NB_i(z)\tag{2}$$

and

$$A(z)=\prod_{i=1}^NA_i(z)\tag{3}$$

So in order to get the coefficients of the polynomials $B(z)$ and $A(z)$ you need to do polynomial multiplication, or, equivalently, convolution of the polynomial coefficients. The coefficients of $B(z)$ are obtained by convolving the coefficients of the individual numerator polynomials $B_i(z)$, and the coefficients of $A(z)$ are obtained in an analogous manner.

The Matlab function sos2tf does exactly that.