We propose subspace methods for three-parameter eigenvalue problems. Such problems arise when separation of variables is applied to separable boundary value problems; a particular example is the Helmholtz equation in ellipsoidal and paraboloidal coordinates. While several subspace methods for two-parameter eigenvalue problems exist, their extensions to a three-parameter setting seem challenging. An inherent difficulty is that, while for two-parameter eigenvalue problems, we can exploit a relation to Sylvester equations to obtain a fast Arnoldi-type method, such a relation does not seem to exist when there are three or more parameters. Instead, we introduce a subspace iteration method with projections onto generalized Krylov subspaces that a...