Time-dependent variation form doesn’t update in the time loop

variational formulation
1

Hey guys. I am calculating how would theta and Vy vary due to time so I construct this code. However, one thing really confuses me is that the solutions are all the same when I try to plot it in 50 different time intervals. I checked a lot but I still don’t understand what is wrong. Any suggestions and help will be appreciated! (I am using fenics2018)

from dolfin import *
import matplotlib.pyplot as plt


# Define 1D mesh and space
mesh = IntervalMesh(10, 0, 1)
Element1 = FiniteElement("CG", mesh.ufl_cell(), 5)
Element2 = FiniteElement("CG", mesh.ufl_cell(), 5)
W_elem = MixedElement([Element1, Element2])
ME = FunctionSpace(mesh, W_elem)


# Define functions
du = TrialFunction(ME)
dtheta, dvy = split(du)
u_old = Function(ME)
u_new = Function(ME)
theta_old, vy_old = split(u_old)
theta_new, vy_new = split(u_new)
vy_mid = 0.5*vy_old + 0.5*vy_new
theta_mid = 0.5*theta_old + 0.5*theta_new
tf = TestFunction(ME)
q, v = split(tf)
f1 = 2 * cos(2 * theta_mid)
f2 = 2 * cos(2 * theta_new)
dt = 1.0


# Boundary condition
def boundary(x, on_boundary):
    return on_boundary


bc = DirichletBC(ME, [Constant(1.57), Constant(0)], boundary)

#Initial Condition
class InitialConditions(UserExpression):
    def eval(self, values, x):
        values[0] = 1.57
        values[1] = 0

    def value_shape(self):
        return (2,)


# Set initial conditions to old function
u_init = InitialConditions()
u_old.interpolate(u_init)

# Weak statement of the equations
F0 = (theta_new-theta_old)/dt*q*dx + theta_mid.dx(0)*q.dx(0)*dx - 0.5*(f1+1)*vy_mid.dx(0)*q*dx
F1 = (-vy_new.dx(0) * v.dx(0)) * dx + f2 * theta_new.dx(0) * v * dx
F = F0 + F1
J = derivative(F, u_new)

#Solver
problem = NonlinearVariationalProblem(F, u_new, bc, J)
solver = NonlinearVariationalSolver(problem)
t = 0.0
T = 50*dt
theta_list = []
vy_list = []

while t < T:
    t += dt
    solver.solve()
    u_old.vector()[:] = u_new.vector()
    theta_list.append(u_new.split()[0])
    vy_list.append(u_new.split()[1])

plot(vy_list[40])
plt.show()