Note that $\cos{\omega t} = (e^{j\omega t} + e^{-j\omega t})/2$ and $\sin{\omega t} = (e^{j\omega t} - e^{-j\omega t})/(2j)$, so your function can be rewritten as $$ f(t) = \sum_{i = 1}^{k}{\frac{a_i - jb_i}{2}e^{j\omega_i t} + \frac{a_i + jb_i}{2}e^{-j\omega_i t}} $$ Obviously, we can solve this function from some differential equation with characteristic roots $\lambda = \pm j\omega_i$ with $i = 1, \cdots, k$. To this end, we can easily construct the characteristic function $$ \Pi_{i = 1}^{k}(\lambda^2 + \omega_i^2) = 0 $$ Expand this equation, replace $\lambda^{i}$ with $f^{(i)}(t)$ and we get the differential equation you want. Initial conditions are determined by values of $a_i$ and $b_i$.