function [x,y,u,it]=laplaciano(f1,f2,g1,g2,n,m,metodo)
% Risolve l'equazione di Laplace sul quadrato di lato 1
% f1, f2 : condzioni sui lati orizzontali
% g1, g2 : condizioni sui lati verticali
% n: discretizzazione sull'asse x
% m: discretizzazione sull'asse y
% metodo: metodo utilizzato per risolvere il sistema lineare
if nargin ==0
    %f1=@(x)zeros(length(x),1);
    f1=@(x)sin(x*pi);
    f2=@(x)sin(x*pi);
    %g1=@(x)zeros(length(x),1);
    g1=@(x)sin(x*pi);
    g2=@(x)sin(x*pi);
end
if nargin <=4
    n=100;
    m=50;
end
if nargin<=5
    metodo='matlab';
end
    
dx=1/n;
dy=1/m;

% MATRICI PIENE
%T=diag(-2*(dx^2+dy^2)*ones(n-1,1))+diag(dy^2*ones(n-2,1),-1)+diag(dy^2*ones(n-2,1),1);
%D=dx^2*eye(n-1);
%A=kron(eye(m-1),T)+kron(diag(ones(m-2,1),-1),D)+kron(diag(ones(m-2,1),1),D);
%As=sparse(A);
%figure(1)
%spy(As)


dim=(m-1)*(n-1); %dimensione del problema
d=-2*(dx^2+dy^2)*ones(dim,1);%diagonale principale
d1=kron(ones(m-1,1),[dy^2*ones(n-2,1);0]);%sottodiagonale
d1(end)=[];% elimino lo zero finale
d2=dx^2*ones(dim-n+1,1);
As= sparse([1:dim]',[1:dim]',d)+sparse([2:dim]',[1:dim-1]',d1,dim,dim)...
   +sparse([n:dim]',[1:dim-n+1]',d2,dim,dim);
As=As+tril(As,-1)';
figure(1)
spy(As)
%termine noto
x=linspace(0,1,n+1)';
y=linspace(0,1,m+1)';
f1x=f1(x(2:n));
f2x=f2(x(2:n));
g1y=g1(y(2:m));
g2y=g2(y(2:m));

b2=[-dy^2*g1y'; zeros(n-3,m-1); -dy^2*g2y'];
b2=b2(:);
b1=zeros((m-1)*(n-1),1);
b1(1:n-1)=-dx^2*f1x;
b1(end-n+2:end)=-dx^2*f2x;
b=b1+b2;


itmax=5000;
tol=1e-12;
it=0;
switch metodo
    case 'matlab'
        u=As\b;
        it=[];
    case 'ilu' %Fattorizzazione LU incompleta
        [L,U]=ilu(As);
        I=eye((m-1)*(n-1));
        B=I-U\(L\As);
        fprintf('La norma della matrice di iterazione è: ||B||=%0.5g \n',norm(B))
        u0=U\(L\b);
        arresto=0;
        it=0;
        % raffinamento iterativo
        while ~arresto && it<itmax
            it=it+1;
            r=b-As*u0;
            d=U\(L\r);
            u1=u0+d;
            arresto=norm(u1-u0)<tol;
            u0=u1;
        end
        u=u1;
    case 'Jacobi'
    case 'Gauss-Seidel'
    case 'rilassamento'
end

u=reshape(u,n-1,m-1)';
u=[[f1(x(1)) f1x' f1(x(n+1))];g1y u g2y; [f2(x(1)) f2x' f2(x(n+1))]];
figure(2)
surf(x,y,u)
shading interp
figure(3)
pcolor(x,y,u)
shading interp
