In solid state physics one encounters a differential equation of the form $$\forall n\in\{1,2...N\}: m\frac{\partial^2 \delta_n(t)}{\partial t²}=\kappa(\delta_{n-1}(t)+\delta_{n+1}(t)-2\delta_n(t))$$ with the boundary condition $\delta_n(t)=\delta_{n+N}(t)$.
What method does one use to find the complete solution to this coupled system of differential equations?
In physics textbooks this equation is often solved by making the Ansatz $\delta_n=\Re\left(a\exp(i(kn-\omega t))\right)$ with $0\le\omega$ by convention. This Ansatz solves the equation iff $$\omega=2\sqrt\frac{\kappa}{m}\left|\sin\left(\frac{k}{2}\right)\right|\quad\textrm{(dispersion relation)}$$ Enforcing the boundary conditions gives a quantization of the parameter $k\phi$ $$\exp(ink)\overset{!}{=}\exp(i(n+N)k)\rightarrow k_j=j\frac{2\pi}{N}, j\in\mathbb{Z}$$ It is not necessary to include all $j$ values though. All distinct solutions of the form of the Ansatz are found if one includes $-N\le j\le N$ (Aliasing). As can be best seen by observing $$\exp(in k_{j})=\exp\left(inj\frac{2\pi}{N}\right)=\exp\left(in(j\textrm{mod}N)\frac{2\pi}{N}\right)\exp\left(inmN\frac{2\pi}{N}\right)=\exp(in k_{j\textrm{mod}N})$$
But how does one gurantee that all solutions are a superpositions of the solutions we found and the solution space is not any larger?
The Ansatz is based on writing $\delta_n(t)$ as it's Fourier series. So as long at that assumption is valid, you have the solution. If the function is continuous, and is bounded (does not diverge to infinity), then the Dirichlet conditions assure you that this is the only solution, since the Fourier transform is unique.