///// Equation de Poisson 1D avec condition Dirichlet //
//// Schema aux differences finies///////
//// Ordre de convergence pour le cas f(x)=sin(%pi x)//
///////////////////////////////////////////////////////
err=[];
for k=2:10
N=2^k;     // valeurs de N: 4, 8, 16,... 1024
dx=1/(N+1);
deff('y=f(x)','y=(%pi^2+1)*sin(%pi*x)')
deff('y=c(x)','y=1')
A=zeros(N,N);
for i=2:N-1
    A(i,i)=2;
    A(i,i+1)=-1;
    A(i,i-1)=-1;
end
A(1,1)=2;A(N,N)=2;A(1,2)=-1;A(N,N-1)=-1;
A=1/(dx^2)*A+eye(N,N);// c=1
//
// aussi: A=toeplitz([2,-1,zeros(1,N-2)]);
b=zeros(N,1);
for i=1:N
    b(i)=f(i*dx);
end
//
u=linsolve(A,-b)
//
deff('y=ue(x)','y=sin(%pi*x)')
x=0:dx:1;uapp=[0;u;0];usol=ue(x);
xex=0:0.001:1;uex=ue(xex);
if (k==4)  then   // trace dans le cas N=4
   figure(1)
   clf()
   plot2d(x,uapp,5)
   plot2d(x,uapp,-1)
   plot2d(xex,uex)
end 
//
err=[err,max(abs(uapp-usol'))];
end
figure(2)  // trace de l'erreur
clf()
plot2d(log(2.^(2:10)),log(abs(err)));
plot2d(log([10^1,10^2]),log([10^(-1),10^(-3)]),2)
legends(['erreur','pente -2'],[1,2])
