Consider the second order linear elliptic differential equation

$$Lu=\sum_{i=1}^d{\partial\over\partial\theta_i}(b_i{\partial\over\partial\theta_i}u)=exp(i\theta_1)$$

This differential equation is defined on the domain $[0,2\pi]^d$ with periodic boundary condition. The coefficients $b_i$ are some polynomial in terms of sines and cosines in $\theta_i$. Is there some algorithm to compute the coefficients of the power series or better yet the Fourier series of the solution $u$?