A methodology is devised for building optimal bases for the generalized Dicke model based on the symmetry adapted variational solution to the problem. At order zero, the matter sector is constructed by distributing Na particles in all the possible two-level subsystems connected with electromagnetic radiation; the next order is obtained when the states of Na − 1 particles are added and distributed again into the two-level subsystems; and so on. In the electromagnetic sector, the order zero for each mode is the direct sum of the Fock spaces, truncated to a value of the corresponding constants of motion of each two-level subsystem; by including contributions of the other modes, the next orders are obtained. As an example of the procedure we consider 4 atoms in the Ξ configuration interacting dipolarly with two modes of electromagnetic radiation. The results may be applied to situations in quantum optics, quantum information, and quantum computing.