I'm trying to understand how to apply the Galerkin method to $u_t - \Delta u = u^3$. I understand how to obtain all of the a-priori estimates using Sobolev embeddings and such but my question concerns the actual discretization procedure where we project onto the finite dimensional subspace spanned by the eigenfunctions of $-\Delta$.

In the linear case we may simply set $u_N = \sum \phi_n c_n$ and plug this in to obtain a set of $N$ O.D.Es which we then show satisfy the same energy bounds. In the non-linear case though we may not just substitute directly because of the $u^3$ term. How can this be dealt with? Is there perhaps a better way to do the approximation?

**Addition:** In this example if we let $u_N = \sum \phi_n c_n^N(t)$ then when we plug this into the weak form of our PDE we obtain and then choose our test function to be one of the basis elements $w_k$ we obtain
$d/dt c_k^N(t) + \sum_{i=1}^n e^{ki}(t) c_k^N(t) = \int (\sum_{n=1}^N c_n^N(t) \phi_n)^3w_k$,
for some coefficients $e^{ki}(t)$.
My question is, how do I deal with the *integral on the right*? I would like to be able to solve this ODE and then say I have a solution.