function [A,B,U]=poiass(gd,m,n,fu)
% Kimmo Kärkkäinen
%
hx=(gd(1,2)-gd(1,1))/(m+1);
hy=(gd(2,2)-gd(2,1))/(n+1);
sxy=hx/hy;syx=hy/hx;

N=m*n;
maindiag=2*(sxy+syx)*ones(N,1);
lowd1=[-syx*ones(m-1,n);zeros(1,n)];
lowd1=lowd(:);

upd1=[0;lowd1(1:N-1)];
fardiag=-sxy*ones(N,1);
A=spdiags([fardiag lowd1 maindiag upd1 fardiag],[(-m),-1 0 1 ,(-m)],N,N);

xk=zeros(n+1,2*(m+1));yk=xk;

....
