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?