pkg load msh
pkg load fem-fenics


fem_init_env ();
problem = 'Evolution';
fem_create_all (problem);

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

V  = Evolution_FunctionSpace (mshd);

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 = Evolution_BilinearForm (V, k);
L = Evolution_LinearForm (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 = Evolution_LinearForm (V, k, u0);
  b = assemble (L, bc);

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

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