function [u]=fdPoisson2d(N,fun,gx0,gy0,gx1,gy1)

% calculates 5-star-FD-approximation of the solution of the Poisson problem on
% [0,1]^2 with (N+1)^2 nodes, uniform mesh width, boundary data g, and rhs
% f; plots discrete solution as surface plot and gives back values of the
% solution at inner nodes

% get x and y values

x=linspace(0,1,N+1);
y=linspace(0,1,N+1);

% define number of nodes

totalnodes=(N+1)^2;
innernodes=(N-1)^2;
outernodes=4*N;

% allocate space

u=zeros(innernodes,1);
g_vec=zeros(outernodes,1);
f_vec=zeros(innernodes,1);
rhs=zeros(innernodes,1);
A=spalloc(innernodes,innernodes,5*innernodes);
B=spalloc(innernodes, outernodes, outernodes-4);

% assemble boundary data

g_vec([1:N+1])=gy0(x);
g_vec(N+2*[1:N-1])=gx0(y(2:end-1));
g_vec(N+1+2*[1:N-1])=gx1(y(2:end-1));
g_vec([end-N:end])=gy1(x);

% assemble right hand side data

f_mat=fun(x([2:N]),y([2:N]));
f_vec=f_mat(:);

% assemble A

E=spdiags(ones(innernodes,2),[-N+1 ,N-1],innernodes,innernodes);
M=spdiags([-ones(N-1,1) 4*ones(N-1,1) -ones(N-1,1)], -1:1, N-1, N-1);

A=N^2*(kron(eye(N-1),M)-E);

% assemble B

I=speye(N-1);
% 4 sides
B([1:N-1],1+[1:N-1])=I;
B(innernodes-N+1+[1:N-1],3*N+[1:N-1])=I;
B([1:N-1]*(N-1),N+1+2*[1:N-1])=I;
B(1+([1:N-1]-1)*(N-1),N+2*[1:N-1])=I;

% assemble right hand side of linear system

rhs=f_vec+N^2*B*g_vec;

% solve linear system

u=mldivide(A,rhs);

% plot solution u

 u_cp=zeros(totalnodes,1);
 % lower boundary
 u_cp(1:N)=g_vec(1:N);
 % upper boundary
 u_cp(N*(N+1)+2:totalnodes)=g_vec(3*N+1:4*N);
 % left boundary
 u_cp(N+2:N+1:totalnodes-N)=g_vec(N+2:2:3*N);
 % right boundary
 u_cp(N+1:N+1:totalnodes-N-1)=g_vec(N+1:2:3*N-1);
 % inner nodes
 u_cp((N+2)*ones(1,(N-1)^2)+kron([0:N-2],(N+1)*ones(1,N-1))+repmat([1:N-1],1,N-1))=u;
 
 % reshape
 u_cp_mat=reshape(u_cp, N+1, N+1);
 
 surf(x,y,u_cp_mat);
 s = sprintf('Surface plot of the discrete solution with mesh width %g',1/N);
 title(s);