Consider the second order linear elliptic differential equation
$$Lu=(\sum_{i=1}^d{\partial^2\over\partial\theta_i^2}+{\partial b_i\over\partial \theta_i}{\partial\over\partial\theta_i})u=exp(i\theta_1)$$$$Lu=(\sum_{i=1}^d{\partial^2\over\partial\theta_i^2}+{\partial b\over\partial \theta_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 coefficientsfunction $b_i$$b$ are some polynomial in terms of sines and cosines in $\theta_k$, $k=1,..,d$. Is there some algorithm to compute the coefficients of the power series or better yet the Fourier series of the solution $u$?