write(%io(2),'Convection-diffusion de polluants ([LuPi], pages 15 et 290)');
//
n=evstr(x_dialog('Entrer le nombre de points interieurs','30'));
nu=evstr(x_dialog('Entrer le coefficient de diffusion','0.06'));
uvent=evstr(x_dialog('vitesse du vent (m/s)','3'));
tfin=evstr(x_dialog('temsp final (s)','10'));
//
h=1/(n+1);npas=tfin/h;k=tfin/npas;u=uvent*ones(n+2,1);
//
//calcul de la matrice a inverser
//Schema de Crank-Nicholson
//condition de Dirichlet
//
v=zeros(1,n-2);
A=toeplitz([(1/k)+(nu/(h^2))-uvent/(2*h),-nu/(2*h^2)+uvent/(2*h),v],[(1/k)+(nu/(h^2))-uvent/(2*h),-nu/(2*h^2),v]');
B=toeplitz([(1/k)-(nu/(h^2))+uvent/(2*h),nu/(2*h^2),v],[(1/k)-(nu/(h^2))+uvent/(2*h),nu/(2*h^2)-uvent/(2*h),v]');
C=inv(A)*B;
//calcul du second membre
deff('y=c(x)','y=x.*(ones(x)-x)');
//
y=c(h:h:n*h)';
ymax=max(y);
//Resolution du systeme
//
for i=1:npas
	y=C*y; xbasc()
      str=string(k*i);
      xstring(0.5,2*ymax/3,['t=',str])
	plot2d((0:h:1)',[0;y;0],1,'011',' ',[0,0,1,0.3]);   
      xtitle('evolution de la concentration','t','C');
end
//