function [x]=laplace(b)

if (size(b,1)==size(b,2) & ~ (size(b,1)==1))
  n=size(b,1);
  b=reshape(b,n^2,1);
else
  n=sqrt(size(b,1))
end;
if (~ (size(b,2)==1 & floor(n)==n))
  error('laplace: Argumentin pitää olla n*n-matriisi tai n^2*1-vektori.');
end;

e=ones(n^2,1);
A=spdiags([-e, -e, 4*e, -e, -e], [-n, -1, 0, 1, n], n^2, n^2);
h=1/(n-1);
b=h*b;
A(1:n,:)=[eye(n) zeros(n,(n-1)*n)];
A(:,1:n)=[eye(n) zeros(n,(n-1)*n)]';
b(1:n)=zeros(n,1);
for i=n:n:(n-1)*n
  A(i:i+1,:)=zeros(2,n^2);
  A(:,i:i+1)=zeros(n^2,2);
  A(i,i)=1;
  A(i+1,i+1)=1;
  b(i)=0;
  b(i+1)=0;
end;
A((n-1)*n+1:n^2,:)=[zeros(n,(n-1)*n) eye(n)];
A(:,(n-1)*n+1:n^2)=[zeros(n,(n-1)*n) eye(n)]';
b((n-1)*n+1:n^2)=zeros(n,1);
x=A\b;
x=reshape(x,n,n);
