// Solve heat equation (22 lines, comments and blank lines included)
// du/dt-Delta u=f 	in Omega
// u=0		  	on the boundary of Omega
// u(t=0)=0

func f=1.;
real T=0.1,dt=0.005,dn=30;

mesh Th=square(dn,dn);

fespace Vh(Th,P1);
Vh u,v,uold;

macro grad(u) [dx(u),dy(u)]//
problem Heat(u,v)=int2d(Th)(u*v/dt+grad(u)'*grad(v))-int2d(Th)(f*v+uold*v/dt)+on(1,2,3,4,u=0);

u=0;	// initialisation
for(real t=0;t<T;t+=dt){
	plot(u,cmm="time="+t);
	Heat;
	uold=u;
}