//
//  equation de Laplace avec conditions de Neuman 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'));
//
alpha=evstr(x_dialog('constante devant u','1'));
//
M=N;
h=1/(N+1);
//
//
//---------------------------------
//  matrice A 
NN=N+2;
L1=-eye(NN,NN);
L2=toeplitz([4 -1 zeros(1,NN-2)]);	//tridiagonale par blocs
//
// corrections avec Neuman
//
L2bis=L2;L2bis(1,2)=-2;L2bis(NN,NN-1)=-2;
//
A=zeros(NN*NN,NN*NN);
for i=1:(NN-1)
  i1=(i-1)*NN;
  i2=i*NN;
  A((i1+1):(i1+NN),(i1+1):(i1+NN)) = L2bis;
  A((i2+1):(i2+NN),(i1+1):(i1+NN)) = L1;
  A((i1+1):(i1+NN),(i2+1):(i2+NN)) = L1;
end;
// corrections avec Neuman
//
A(1:NN,1:NN)=L2bis;
A(1:NN,NN+1:2*NN)=2*L1;
//
i1=(NN-1)*NN;
i2=i1-NN;
A((i1+1):(i1+NN),(i1+1):(i1+NN)) =L2bis;
A((i1+1):(i1+NN),(i2+1):(i2+NN)) =2*L1;
//
//
A=1/(h^2)*A+alpha*eye(NN^2,NN^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(NN*NN,1);
for j=1:NN
 for i=1:NN
      b((j-1)*NN+i)=f((i-1)*h,(j-1)*h)  
  end
end
//----------------------------------
// resolution de AU=b
//
U=inv(A)*b;
//----------------------------------
// reecriture sous forme matrice de U (avec bords nuls)
//
umat=matrix(U,NN,NN); // reecriture colonne par colonne (adaptée à plot3d)
//
//----------------------------------
//  tracé
clf()
plot3d(0:h:1,0:h:1,umat)       
