
% MODIFIED VERSION: MARCH 13, 2008.

demo_example=3;

switch demo_example
case 1
    f=inline('zeros(size(x))','x','y');
    g_down=inline('ones(size(x))','x','y');
    g_up=inline('ones(size(x))','x','y');
    g_left=inline('ones(size(x))','x','y');
    g_right=inline('ones(size(x))','x','y');    
    
    solution=inline('ones(size(x))','x','y');
    
case 2
    f=inline('zeros(size(x))','x','y');
    g_down=inline('exp(pi*x)','x','y');
    g_up=inline('-exp(pi*x)','x','y');
    g_left=inline('cos(pi*y)','x','y');
    g_right=inline('((exp(1))^pi)*cos(pi*y)','x','y');    
    
    solution=inline('(exp(pi*x)).*cos(pi*y)','x','y');
    
case 3
    f=inline('(-2*(pi^2))*sin(pi*x).*sin(pi*y)','x','y');
    g_down=inline('zeros(size(x))','x','y');
    g_up=inline('zeros(size(x))','x','y');
    g_left=inline('zeros(size(x))','x','y');
    g_right=inline('zeros(size(x))','x','y');
    
    solution=inline('(sin(pi*x)).*sin(pi*y)','x','y');
    
otherwise
    f=inline('ones(size(x))','x','y');
    g_down=inline('zeros(size(x))','x','y');
    g_up=inline('zeros(size(x))','x','y');
    g_left=inline('zeros(size(x))','x','y');
    g_right=inline('zeros(size(x))','x','y');
    
    solution=inline('(sin(pi*x)).*sin(pi*y)','x','y');
end

for index=2:5
    n=2^index;
    
    h=1/(n+1); x=(h:h:1-h)'; y=x;
    [X,Y]=meshgrid(x,y);
    X=X'; Y=Y';
    
    u=poisson5pts(n,f,g_left,g_right,g_down,g_up);
    
    % USO RESHAPE COSI' LA SOLUZIONE HA LE STESSE DIMENSIONI DELLE MATRICI X, Y.
    Z=(reshape(u,n,n));
    
    if demo_example <=3
        
        V=feval(solution,X,Y);
        
        err(index)=norm(V(:)-Z(:),inf);
        if index == 1
            fprintf('\n \t [n]: %4.0f [ERR]: %2.2e',n,err(index));
        else
            fprintf('\n \t [n]: %4.0f [ERR]: %2.2e [RATIO]: %2.2f',n,err(index),err(index-1)/err(index));
        end
        
    end
end

surf(X,Y,Z);