%
% Poissonin yhtälö yksikköneliössä (säännöllinen hila)
%
function poisson_dir4(nn)
global n;
n=nn;
h=1/(n+1);
h2=h*h;
xx=0:h:1;
yy=0:h:1;
[X,Y]=meshgrid(xx,yy);
f=X.*Y;
u=zeros(n+2,n+2);
b=zeros((n+2)*(n+2),1); % oikean puolen vektori
for j=1:n+2
    for i=1:n+2
        ii = (j-1)*(n+2)+i;
        if i==1 | i==n+2 | j==1 | j==n+2
            %
            %   Jos hilapiste on reunalla, niin sitä
            %   vastaa yhtälö u(ii)=0.0
            %
            b(ii)=0.0;
        else
            %
            % Muussa tapauksessa yhtälö on muotoa
            %  4*u(ii)-u(ii-1)-u(ii+1)-u(ii+n+2)-uu(ii-n-2)=h*h*f(j,i)
            %
            b(ii)=h2*f(j,i);
        end
    end
end
%
% Ratkaistaan liittogradienttimenetelmällä yhtälöryhmä A*uu=b:
%
uu = pcg(@matvec,b,1e-12,200);
%
% Siirretään ratkaisun hilapistearvot vektorista uu matriisiin u:
%
for j=1:n+2
    for i=1:n+2
        u(j,i)=uu((j-1)*(n+2)+i);
    end
end
max(max(u))
contourf(X,Y,u,10); axis equal;
colorbar('vert')

%
% Funktio palauttaa matriisi-vektoritulon A*v arvon annetulle v.
%
function av=matvec(v)
global n;
av=zeros((n+2)*(n+2),1);    %  matvec should return column vector
for j=1:n+2
    for i=1:n+2
        ii = (j-1)*(n+2)+i;
        if i==1 | i==n+2 | j==1 | j==n+2
            %
            %   Jos hilapiste on reunalla, niin sitä
            %   vastaa yhtälö u(ii)=0.0
            %
            av(ii)=v(ii);
        else
            %
            % Muussa tapauksessa yhtälö on muotoa
            %  4*u(ii)-u(ii-1)-u(ii+1)-u(ii+n+2)-uu(ii-n-2)=h*h*f(j,i)
            %
            av(ii)=4*v(ii)-v(ii-1)-v(ii+1)-v(ii+n+2)-v(ii-n-2);
        end
    end
end