from dolfin import *

# Mesh and function space
mesh = UnitSquare(8, 8)
V = FunctionSpace(mesh, "Lagrange", 1)

# Time variables
dt = Constant(0.3); t = float(dt); T = 1.8
g_expr = '1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t'
g = Expression(g_expr, alpha=3.0, beta=1.2, t=0,
               degree=2)

# Previous and current solution
u0 = interpolate(g, V)
u1 = Function(V)

# Variational problem at each time
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(1.2 - 2. - 2*3.0)
a = u*v*dx + dt*inner(grad(u), grad(v))*dx
L = u0*v*dx + dt*f*v*dx
bc = DirichletBC(V, g, "on_boundary")

while (t <= T):
    # Solve
    g.t = t
    solve(a == L, u1, bc)

    # Update
    u0.assign(u1)
    t += float(dt)

plot(u1, title="Approximated final solution")
interactive()