// Equation de la chaleur - schema implicte
clear; close,
da=gda(); da.auto_clear='on';

V=1;	// vitesse de convection
T=0.1;	// temps final
cfl=10.;
	
N=[20,40,80,160, 320, 640];	// Nb de point de discretisation

erreur=[];
for n=N
	dx=1/n;
	dt=cfl*dx^2/2.
//	dt=dx;
	x=linspace(0,1,n+1);
	//u0=(1.-4.*((x-0.5)^2))';
	u0=cos(%pi*x)';
	u=u0;
	plot(x,u);
	A=diag(sparse(2.*ones(1,n+1)))-diag(sparse(ones(1,n)),1)-diag(sparse(ones(1,n)),-1);
	A(1,2)=-2.;
	A(n+1,n)=-2.;
	A=A/(dx)^2;
	B=speye(n+1,n+1)+dt*A;
	h=lufact(B);
	for t=dt:dt:T,
		u=lusolve(h,u);//schéma implicte
		//	plot2d(x,[u u0]);
	end;
	uexact=exp(-%pi^2*T)*cos(%pi*x)';
	erreur=[erreur,sqrt(((u-uexact)'*(u-uexact)*dx))];
//	plot2d(x,[u u0 uexact]);

end;
dN=ones(N)./N;
dN2=ones(N)./(N^2);
plot2d(N,[erreur',dN',dN2'],logflag='ll');