Error in computing solve() when we are using u at previous time steps

variational formulation
1

Hello everyone

Here I want to compute value of solution u using data u at t =0.2 so, first I append all the solution in a list U and try to compute u at t=0.6 using u at t=0.2 in second for loop but it shows error.
Here u_n is the solution u at t=0 and f is given function.
Please suggest something in order to correct it.

from fenics import *
import numpy as np

T = 1.0 # final time
num_steps = 5 # number of time steps
dt = T / num_steps # time step size
alpha = 3 # parameter alpha
beta = 1.2 # parameter beta
s = 0.4 # delay
k = s / dt


# Create mesh and define function space
nx = ny = 8
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition
u_D = Expression('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t', degree=2, alpha=alpha, beta=beta, t=t)
def boundary(x, on_boundary):
    return on_boundary
bc = DirichletBC(V, u_D, boundary)

# Define initial value
u_n = interpolate(u_D, V)
#u_n = project(u_D, V)
# Define variational proble
u = TrialFunction(V)
v = TestFunction(V)
f = Expression('beta - 3 -2*alpha - x[0]*x[0] - alpha*x[1]*x[1] - beta*(t-s)', degree=2, alpha=alpha, beta=beta, s=s, t=0)
phi = Expression('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t', degree=2, alpha=alpha, beta=beta, t=-s)

t=0
U =[u_n]
for i in range(num_steps):
    F1 = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx - dt*phi*v*dx 
    a1, L1 = lhs(F1), rhs(F1)
    u = Function(V)
    if i <= k:
       #Update current time
       t = (i+1)*dt
       u_D.t = t
       f.t = t
       phi.t = t
       #Compute solution
       solve(a1 == L1, u, bc)
       U.append(u)
       #print(u_nodes.append(u))
       u_n.assign(u)
       break

print(U)
for i in range(num_steps):
    F2 = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx - dt*U[i-k]*v*dx 
    a2, L2 = lhs(F2), rhs(F2)  
    u = Function(V)
    if i > k:
       t = (i+1)*dt
       u_D.t = t
       f.t = t
       phi.t = t
       #Compute solution
       solve(a2 == L2, u, bc)
2

See: Simple heat problem also doesn't vary due to time - #2 by nate