To answer my question for future reference, I think the issue is that dPi is the integral of weak form of the Euler-Lagrange equation of Pi, while u0 is a Function of space.
So, I need to use u0_int = dot(u0, v)*dx and then use u1 = -deltax * dPi_val + u0_int
2 Likes