clear
L=1;

J=30;
h=L/(J+1);
k=h^2/2;
N=floor(0.15/k)//nb de pas de temps

a=0.4;// pour la CI
deff('y=u0(x)','y=bool2s((x<=a))');

Ah=2*diag(ones(J,1),0)-diag(ones(J-1,1),1)-diag(ones(J-1,1),-1);
Ah=Ah/h^2;
Fh=zeros(J,1);Fh(1)=1/h^2;

x=[h:h:L-h];
U=u0(x');
figure(0)
clf
  plot2d([0 x 1],[1 U' 0])
    halt()
TabU=zeros(J+2,N+1);
TabU(1,1:N+1)=ones(1,N+1);
t=0;
tabt=t;
for n=1:N
    t=t+k;
    tabt=[tabt t];
    U=U-k*Ah*U+k*Fh;
    TabU(2:J+1,n+1)=U;
    //FIGURE
    plot2d([0 x 1],[1 U' 0])
    xtitle('','x','u')
    halt()
end
figure(1)
plot3d([0 x 1],tabt,TabU)

