We present the analytical solution of the Tavis-Cummings (TC) model for more than one qubit inhomogeneously coupled to a single mode radiation field beyond the rotating-wave approximation (RWA). The significant advantage of the displaced oscillator basis enables us to apply the same truncation techniques adopted in the single qubit Jaynes-Cummings (JC) model to the multiple qubits system. The derived analytical spectrum match perfectly the exact diagonalization numerical solutions of the inhomogeneous TC model in the parameter regime where the qubits transition frequencies are far off-resonance with the field frequency and the interaction strengths reach the ultra-strong coupling regime. The two-qubit TC model is quasi-exactly solvable because part of the spectra can be determined exactly in the homogeneous coupling case with two identical qubits or with symmetric(asymmetric) detuning. By means of the fidelity of quantum states we identify several nontrivial level crossing points in the same parity subspace, which implies that homogeneous coupled two-qubit TC model with $omega_1=omega_2$ or $omega_1pmomega_2=2omega_c$ is integrable. We further explore the time evolution of the qubits population inversion and the entanglement behavior taking two qubits as an example. The analytical methods provide unexpectedly accurate results in describing the dynamics of the qubit in the present experimentally accessible coupling regime, showing that the collapse-revival phenomena emerge, survive, and are finally destroyed when the coupling strength increases beyond the ultra-strong coupling regime. The suggested procedure applies readily to the multiple qubits system such as the GHZ state entanglement evolution and quantum entanglement between a single photon and superconducting qubits of particular experiment interest.