//
//  equation de Laplace avec conditions de Dirichlet sur [0,1]^2
// 
//  N.B. les noeuds sont numérotés ligne par ligne à partir du bas gauche
//
N=evstr(x_dialog('nombre de points intérieurs','10'));
fval=evstr(x_dialog('valeur de f(x) au centre du domaine','10'));
//
c=evstr(x_dialog('constante devant u','1'));
M=N;
h=1/(N+1);
//---------------------------------
//  matrice A (cas M=N)
L1=-eye(N,N);
L2=toeplitz([4 -1 zeros(1,N-2)]);	//tridiagonale par blocs
A=zeros(N*N,N*N);
for i=1:(N-1)
  i1=(i-1)*N;
  i2=i*N;
  A((i1+1):(i1+N),(i1+1):(i1+N)) = L2;
  A((i2+1):(i2+N),(i1+1):(i1+N)) = L1;
  A((i1+1):(i1+N),(i2+1):(i2+N)) = L1;
end;
i1=(N-1)*N;
A((i1+1):(i1+N),(i1+1):(i1+N)) =L2;
A=1/(h^2)*A;
//
A=A+c*eye(N^2,N^2);
//-------------------------------
// second membre
//
function fu=f(x,y)
    fu=0;
    if (x>1/3)&(x<2/3)&(y>3/5)&(y<4/5) then
        fu=fval;
    end    
endfunction
b=zeros(N*N,1);
for j=1:N
 for i=1:N
      b((j-1)*N+i)=f(i*h,j*h)  
    end
end
//----------------------------------
// resolution de AU=b
//
U=inv(A)*b;
//----------------------------------
//
// reecriture sous forme matrice de U (rajout des bords nuls)
//
umat=zeros(N+2,N+2)
for j=2:N+1
    for i=2:N+1
        umat(i,j)=U((j-2)*N+i-1)  // colonne par colonne (adaptée à plot3d)
    end
end
//----------------------------------
//  tracé
clf()
plot3d(0:h:1,0:h:1,umat)       
