Consider \(m\)-fold products \(\eta = \eta_1 \dots \eta_m\) of special functions. Each of the factors satisfies a differential equation (of order \(\geq 1\) \([\geq 2\), in general]). A general method is introduced so that this product is the first component of a solution of a first order system of differential equations (Theorem 2.1). The construction of the system is explicitly given, but too lengthy to be reproduced here. If each of the functions \(\eta_j\) is a solution of an eigenvalue problem (with common eigenvalue parameter for all factors), then the system for the products can also be written as an eigenvalue problem, and the question of eigenfunction expansions for the product arises. A particular example of a product of Bessel functions \(J_{\nu_1 + {k \over 2}} (az)\) and Whittaker functions \(M_{\kappa, \nu_2 + {k \over 2}} (bz)\) is considered, where \(k \in \mathbb{Z}\). It is shown that holomorphic functions \(f\) on the Riemann surface of the logarithm satisfying the Floquet condition \(f(ze^{2 \pi i}) = e^{\nu_1 + \nu_2 + {1 \over 2}} f(z)\) can be expanded into series of \(J_{\nu_r + {k \over 2}} (az) M_{\kappa, \nu_s+ {k\over 2}} (bz)\), \((r,s) \in \{(1,2), (2,1)\}\), \(k \in \mathbb{Z}\). It is expected that this will become a standard method to obtain expansions into products of special functions.