Why am I getting only the trivial solution to pair of nonlinear PDEs?

[dokken]

@Joshua_Goldberg

u_n = interpolate(u_0, V)
                 
# Define variational problem

u = Function(V) # current step
u_n = Function(V) # previous step

You are redefining u_n, so as a fist step, i would suggest removing the last line of the code i quoted above, such that you only have

u_n = interpolate(u_0, V)
u_n1, u_n2 = split(u_n)