% Harj. 7 teht. 3, Poissonin yhtälö neliöalueessa, 0-reunaehdot
function [u,X,Y]=Poisson0(nelio,f,n)
% 
% Syötteet: nelio   --  [a,b,s] : vasen alanurkka, sivun pituus
% "Lähdefunktio f"  --  2:n muuttujan funktio
%           n       --  sisäjakopisteiden lkm. (yhdessä dim:ssa)


a=nelio(1);b=nelio(2); s=nelio(3);
if nargin < 3, n = 31; end
h = s/(n+1);
x=a:h:a+s; y=b:h:b+s;
[X,Y] = meshgrid(x,y);
F=feval(f,X,Y);
e = ones(n,1);
D = spdiags([e -2*e e],[-1 0 1],n,n);
I = eye(n,n);
A = (1/h^2)*(kron(D,I) + kron(I,D));
sisa_i = 2:n+1;
F0 = F(sisa_i,sisa_i);
U = A\F0(:);
u = zeros(n+2,n+2);
u(sisa_i,sisa_i) = reshape(U,n,n);
