function [A,b]=potenziale(nx,ny,rho);
%
% Costruzione della matrice A e del termine noto b per l'approssimazione 
% del potenziale elettrico su un quadrato
%
% Input: nx, ny= numero di nodi (interni al quadrato) di discretizzazione
%                lungo la direzione x e y, rispettivamente
%        rho =  macro contenente la funzione densita' di carica,
%               definita come rho=@(x,y)[...] su variabili array
%
% Output: A= matrice in formato sparse, di dimensione n=nx*ny
%         b= vettore (colonna) termine noto, di dimensione n
%
hx=1/(nx+1);hy=1/(ny+1);
x1=linspace(hx,1-hx,nx);
y1=linspace(hy,1-hy,ny);
A=lapdf(nx,ny,hx,hy);
b=1.128823459953265e+11*valutaf(x1,y1,rho);

function [A]=lapdf(nx,ny,hx,hy);
%
% costruzione di matrice per il laplaciano, 
% differenze finite centrate, con matrice sparsa.
%
n=nx*ny;
hx2=1/(hx*hx);hy2=1/(hy*hy);
alpha=2*(hx2+hy2); beta=-hx2; gamma=-hy2;

% costruzione di A

e=ones(n,1);

A=spdiags([beta*e alpha*e beta*e], -1:1, n,n);
A=A+spdiags([gamma*e gamma*e], -nx:2*nx:nx, n,n);

% devo toglie gli elementi dei punti vicini ai vertici.


for j=1:ny-1; k=j*nx; A(k,k+1)=0; end
for j=2:ny; k=(j-1)*nx+1; A(k,k-1)=0; end


function [b]=valutaf(xx,yy,f);
%
% Valutazione di una funzione matematica nei nodi (x,y)
%
nx=length(xx); ny=length(yy);
n=nx*ny; b=zeros(n,1);
for j=1:ny
for i=1:nx
k=(j-1)*nx+i; b(k)=f(xx(i),yy(j));
end
end


