pkg load fem-fenics msh

problem = 'Evolution';
import_ufl_Problem (problem);

mesh = Mesh (msh2m_structured_mesh (0:0.05:1, 0:0.05:1, 1, 1:4));

V = FunctionSpace(problem, mesh);

bc = DirichletBC (V, @(x,y) 1, 1:4);

t = 0;
dt = 0.01;
T = 0.3;

k = Expression ('dt', @(x,y) dt);

u0 = Expression ('u0', @(x,y) 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02));

a = BilinearForm (problem, V, V, k);
L = LinearForm (problem, V, k, u0);

A = assemble (a, bc);

# solve the problem for each time step
# We need to update only the lhs
while t < T
  t += dt;

  # we can pass u0 to the lhs indifferently as a fem_coeff or
  # as a fem_func
  L = LinearForm (problem, V, k, u0);
  b = assemble (L, bc);

  u = A \ b;
 
  u0 = Function ('u0', V, u);

  #press Q to make the plot continue
  plot (u0);
end