function u = poissonSolver2D(f);

[M,N] = size(f);

Tm = sparse(diag(-4*ones(M,1))+diag([ones(M-1,1)],1)+diag([ones(M-1,1)],-1));

lambdam = diag(idst(eye(M))*Tm*dst(eye(M)));

c = idst(eye(M))*f;

v = zeros(size(f));

disp('Starting solution.');

for r = 1:M
    if ~mod(r,ceil(M/10))
        disp(['Computation ' num2str(floor(100*r/M)) '% complete.']);    
    end
    
    v(r,:) = [tridisolve(ones(N-1,1),repmat(lambdam(r),N,1),ones(N-1,1),[c(r,:)]')]';
end

u = dst(eye(M))*v;
