This paper is in a certain sense a continuation of a previous one by the same authors [J. Math. Anal. Appl. 179, No. 2, 610-629 (1993; Zbl 0802.30001)]. Instead of axial functions here bi-axial functions \(f_{k,1}\) are considered. These are solutions of the equation \(\partial_{\vec x} f_{k,1} = 0\), where \(\partial_{\vec x}\) is the Dirac operator of Euclidean \(m\)-space, of the form \[ f_{k,1} (x) = \bigl [A_{k,1} (\rho_ 1, \rho_ 2) + \vec \omega_ 1 B_{k,1} (\rho_ 1, \rho_ 2) + \vec \omega_ 2 C_{k,1} (\rho_ 1, \rho_ 2) + \vec \omega_ 1 \vec \omega_ 2 D_{k,1} (\rho_ 1, \rho_ 2) \bigr] P_{k,1} (\vec \omega_ 1, \vec \omega_ 2) \] where \(\vec x = \vec x_ 1 + \vec x_ 2\), \(\vec x_ i = \rho_ i \vec \omega i\), \(i = 1,2\), \(\vec \omega_ 1\) is a unit vector in \(\mathbb{R}^ p\), \(\vec \omega_ 2\) in \(\mathbb{R}^ q\) \((p + q = m)\), \(P_{k,1} (\vec \omega_ 1, \vec \omega_ 2)\) is a spherical (i.e. homogeneous) monogenic in both its first and its second variable, and \(a_{k,1}, \dots\) are scalar. This generalises the case \(q = 1\) in the previous paper. The central result is that, given a holomorphic function \(f\), then \[ \begin{multlined} \int_{-1}^ 1 \int_{-1}^ 1(1 - t^ 2)^{k + (p - 3)/2} (1 - s^ 2)^{1 + (q - 3)/2} f(\rho_ 1t + i \rho_ 2 s) \\ (1 \vec \omega_ 1t + i \vec \omega_ 2 s + i \vec \omega_ 1 \vec \omega_ 2 st) dt ds P_{k,1} (\vec x_ 1, \vec x_ 2) \end{multlined} \] is a bi-axial monogenic function, where the integration path can be deformed in order to avoid singularities of \(f\). This construction is compared with a construction using plane waves (monogenic functions depending only on the projection of \(x\) on a fixed plane). Various applications are given. Mainly there is a study of the case where \(f(z) = z^ \alpha\), first for integer \(\alpha\), later for the general case. Further the axial exponential function, where \(f(z) = ce^ z\) is treated \((c\) is a normalising constant), where an estimate for the modulus of the function is given. A second topic of the paper is the theory of generalised orthogonal (Hermite and Gegenbauer) polynomials. These polynomials arise in the series expansion of bi-axial monogenic functions. There is a difference however between the case \(q = 1\) and the general case. If \(q = 1\) it is sufficient to know the function for \(\vec x_ 2 = 0\); because of the extra degrees of freedom for \(p > 1\), it is necessary to describe the asymptotic behaviour of the function as \(\rho_ 2 \to 0\). Starting from the standard generating function both Hermite and Gegenbauer polynomials are calculated, and it is proved that they are (up to normalising factors) the same as the ones described in the previous paper.