We are dealing with Poisson's equation with inhomogeneous BCs and (according to the text), any geometry. $$\nabla^2 u = Q(x,y); 0 \le x \le L,0\le y \le H \\ u = \alpha \mid x,y\in\partial\Omega$$
The full derivation for the the Fourier coefficient is $b_i = -\frac{1}{\lambda_i}\frac{\int\int\phi_iQ + \oint \alpha \nabla \phi_i \cdot \hat{n} dxdy}{\int\int \phi_i^2dxdy}$ where $u = \alpha$ on the boundary. The book I'm reading explains how to arrive to this well enough, but I'm confused on how to get to one of the equations at an earlier step in the derivation, which is:
$$\tag{1} \boxed{b_i = \frac{\int\int u\phi_idxdy}{\int\int \phi_i^2 dxdy}}$$
To preface this, we have that the solution is inhomogeneous on the boundaries, and we choose to use $\phi$ such that it is homogeneous on the boundary (My guess is that we use $\phi = \sin(\frac{n\pi x}{L})\sin(\frac{m\pi y}{H})$ since it is homogeneous on the boundaries). From this, we expand the eigenfunction to define u as:
$$\tag{2} u = \sum_i b_i \phi_i$$
and because $\phi = 0 \mid x,y\in\partial\Omega$, it is no longer true that $\nabla^2 u = \sum_i b_i \nabla^2\phi_i$ on the boundary, since on the boundary we would get: $$\nabla^2 \alpha \ne \sum_i b_i\nabla^2 0$$
Immediately from these givens, the book then says that we can use (2) to get (1). I tried deriving this myself, but I consistently get $Q(x,y)$ in the integrand, not $u$.
Just one additional bit of context, it still holds that: $$\tag{3} \nabla^2 \phi_i = -\lambda_i\phi_i$$
Here is my process:
plugging (2) into the PDE:
$$\tag{4} \nabla^2 u = Q$$ gives
$$\tag{5} \sum_i b_i \nabla^2 \phi_i = Q$$ and maintaining that (3) is true gives:
$$\tag{6} -\sum_i\lambda_i b_i \phi_i = Q$$ and now we're at the step where we apply the
orthogonality of sines to isolate $b_i$. From the orthogonality of sines:
$$\tag{7} \int_0^L \sin(\frac{n\pi x}{L}) \sin(\frac{p\pi x}{L}) = \cases{0 & n=/= p \\ L/2 & n = p}$$
and substituting $\sin(n\pi x/L)$ with $\phi_i$ and $\sin(p\pi x/L)$ with $\phi_j$, and lastly understanding that $\int a dx \int b dy = \int\int ab dxdy$ lets us easily extend this property to 2-D and say (presumably) that $\phi_i = \sin(\frac{n\pi x}{L})\sin(\frac{m \pi y}{H})$ and $\phi_j = \sin(\frac{p\pi x}{L})\sin(\frac{q \pi y}{H})$ and this would yield $\frac{LH}{4}$ at $n=p,m=q$ and 0 otherwise. We can multiply both sides of (6) by $\phi_j$ and integrate to get:
$$\tag{8} -\sum_i \lambda_i b_i \int_0^H\int_0^L \phi_i \phi_j dxdy = \int_0^H\int_0^L Q \phi_j dxdy$$
And applying (7) (not applying the result, using this to eliminate the $\sum_i$) we get:
$$\tag{9} -\lambda_jb_j\int_0^H\int_0^L \phi_j^2 dxdy= \int_0^H\int_0^L Q\phi_j dxdy $$
And lastly, we just divide the integrals and $\lambda_j$ over, and arbitrarily removing the bounds of integration and changing j to i to match our notation to that of (1), to get:
$$\tag{10} b_i = -\frac{1}{\lambda_i}\frac{\int\int Q\phi_idxdy}{\int \int \phi_i^2 dxdy}$$
So just from looking at (10) and and how it differs from (1), I'm guessing $$\tag{11} -\frac{Q}{\lambda_i} = u $$, but I don't see how to arrive to this conclusion.