Consider the second order linear elliptic differential equation

$$Lu=(\sum_{i=1}^d{\partial^2\over\partial\theta_i^2}+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$?