//
// Resolution de l'equation 
//  du/dt -Laplacien(u)=f.
// schema de Gear


mesh Th=square(30,30,[-1+2.*x,-1+2.*y]);

func f=1.;

real dt=0.1;
real T=2; //temps final

fespace Vh(Th,P1);

Vh u,u1,u0,v;

problem chaleur(u,v,solver=LU)=
	int2d(Th)(
	(dx(u)*dx(v)+dy(u)*dy(v))
	+3*u*v/dt
	)
	-int2d(Th)(
	f*v
	+(4.*u1-u0)*v/dt
	)
	+on(1,2,3,4,u=0)
;



u0=0.;
u1=0.;
string legende="t=0";
plot(u,fill=1,cmm=legende,wait=1);

int Niter=T/dt;

for(int i=1;i<Niter+1;i++){
chaleur;
u0=u1;
u1=u;
string legende="t="+dt*i;
plot(u,fill=1,cmm=legende);
}
legende="temps tinal, t="+dt*Niter;
plot(u,fill=1,wait=1,cmm=legende);