Problem: codying a variational formulation

Hello, I am a newbie on FreeFEM and I want to solve a problem. I have this variational formula (written in LaTeX code) :
\begin{align*}
\int_\Omega \frac{du}{dt}v \text{ } dx + \int_\Omega k(x)\nabla u \cdot\nabla v \text{ } dx = \int_\Omega \rho u v(1-u^2) \text{ } dx.
\end{align*}.

I have written in FreeFEM this code:
real rho = 1.0;
real k = 1.0;
real dt = 0.1;
func u0 = 0.5; // Condición inicial

mesh Th = square(20, 20); // Domínio Omega, cuadrícula 20x20
fespace Vh(Th, P1); // Espacio de funciones P1

Vh u=u0, v, uold; // Funciones de prueba, de base y solución inicial

// Condiciones de contorno

real T = 10; // Tiempo final

// Formulación variacional
problem Ejercicio(u, v,solver=LU)
= int2d(Th)( u/dtv) // Derivada temporal
+ int2d(Th)(k
( dx(u) * dx(v) + dy(u) * dy(v)))// Término de difusión
- int2d(rhouv*(1 -u*u)) // Término no lineal
+ on(1, u=0); // Condición de Dirichlet en el borde

for(real t = 0; t < T; t += dt){
Ejercicio;// uold ->u //here the elliptic problem is solved
uold = u; // update; equivalent to u^{n-1} = u^n
plot(u);
}

but it doesn´t work because in the lineal term ( - int2d(rhouv*(1 -uu))) FreeFEM doesn’t know how treat the product uu (it interprets it as a matrix product). I’m not sure if I need to use another solver or some other way. Thank you very much in advance.