// Equation de la chaleur
//clear; close,
//da=gda(); da.auto_clear='on';
// instabilité schéma explicite
//T=0.1;	// temps final
//cfl=1.01;
T=0.01
cfl=1.00
	theta=0.
//N=[20,40,80,160 ]';	// Nb de point de discretisation
J=input("nombre de points en espace ")

	dx=1/J;
	dt=cfl*dx^2/2.
disp(dt,"dt")
disp(T/dt,"iter")
N=floor(T/dt)
disp(N,'nombre d iterations')
Nt=input("nombre d instantanés pour la visu ")
pasN=N/(Nt-1)
disp(pasN,"pas en temps des visus de sortie ")
// W_d=spzeros(J+1,Nt)

count=0
pascount=0
nn=1
	x=linspace(0,1,J+1);
function y=u0(x0)
// données régulières
//    y=(1.-4.*((x0-0.5).^2))
// y=cos(J*%pi*x0)
	y=cos(5*%pi*x0)';
// donnees chocs
// y=double(x1<0.25)+double(x1>0.75)
endfunction
w=feval(x,u0)
	u=w';

plot(x,u');
	A=diag(sparse(2.*ones(1,J+1)))-diag(sparse(ones(1,J)),1)-diag(sparse(ones(1,J)),-1);
	A(1,2)=-2.;
	A(J+1,J)=-2.;
	A=A/(dx)^2;
B=speye(J+1,J+1)+theta*dt*A;

	B2=speye(J+1,J+1)-(1.-theta)*dt*A;
for t=dt:dt:T,
       count=count+1
		v=B2*u;
		u=B\v;// theta-schema	
//v=u-dt*A*u;//schéma explicite
//			v=B\u;//schéma implicte
	//u=v;
    if (count==floor(pasN)) then
        pascount=pascount+1
        disp(pascount*floor(pasN),'iteration')
        count=0
        nn=nn+1
       
      plot2d(x,u);
    end

	end;
scf
// si pas solution excate 
//plot2d(x,u)		
// uniquement pour une donnée sin ou cos
	uexact=exp(-25*%pi^2*T)*cos(5*%pi*x)';
	plot2d(x,[u   uexact]);


