function FpointNeu1 = FpointNeu1(n,plot)

% calculating the numerical solution of

%         -laplacian operator of u = f in [0,1]x[0,1]

% 5point-star
% von Neumann boundary conditions
% cartesian grid with element length 1/n
% lexicographic numbering

% f =4pi sin(2piX)( pi cos(2pi y^2)(1+4y^2)+sin(2pi y^2) )
% => the analytic solution u0 is:

%        u0=sin(2piX)cos(2piY^2)


%************* the von Neumann boundary conditon t ***************

% the values of g are deduced of the known analytic solution u0
x=1/n*(1:n-1);

% lower boundary to: y=0  du/dn=-du/dy
% -du0/dy =-4piY sin(2pi X)sin(2pi Y^2)
to=zeros(1,n-1);

% left boundary tl: x=0  du0/dn=-du0/dx
%-du0/dx =-2pi cos(2pi X)cos(2pi Y^2)
tl=-2*pi*cos(2*pi*x.^2);

% right boundary tr: x=1 du0/dn=du0/dx
% du0/dx =2pi cos(2pi X)cos(2pi Y^2)
tr=-tl;

% upper boundary tu: y=1 du/dn=-du/dy
% du0/dy =-4piY sin(2pi X)sin(2pi Y^2)
tu=zeros(1,n-1);

%**************  vectors f and p ****************

f=[];
for j=1:n-1
   f=[f,sin(2*pi*x)*4*pi*(pi*cos(2*pi*(j/n)^2)+sin(2*pi*(j/n)^2)+4*pi*(j/n)^2*cos(2*pi*(j/n)^2))];
end

% p will be the right side of the finite difference equation A*u=p
p=f;

% when the star is lieing next to the border, we have to change the rigth side p

% lower border
p(1:n-1)                =p(1:n-1)                +n*to;
% left border
p(1:n-1:(n-1)*(n-2)+1)  =p(1:n-1:(n-1)*(n-2)+1)  +n*tl;
% right border
p(n-1:n-1:(n-1)^2)      =p(n-1:n-1:(n-1)^2)      +n*tr;
% upper border
p((n-1)*(n-2)+1:(n-1)^2)=p((n-1)*(n-2)+1:(n-1)^2)+n*tu;

p=p';

% assure that p is in the image of A

v=ones((n-1)^2,1); %is in the kernel of A
p=p-(v'*p)/(v'*v)*v; 

% Dirichlet line of A
p(1)=sin(2*pi/n)*cos(2*pi/n^2);  %=> u(1)

%***************  Matrix A **********************

e=ones(n-1,1);

S=spdiags([e -2*e e],[-1 0 1],n-1,n-1);
S(1,1)=-1;
S(n-1,n-1)=-1;

I=spdiags([-e],[0],n-1,n-1);

A=kron(I,S)+kron(S,I);
A=A*n^2;

% Dirichlet-line
A(1,1)=1;
A(1,2)=0;
A(1,n)=0;

%*************** solution u **************************

u=A\p;
u=u';

%****************  comparing solutions  ****************

% the analytic solution as a vector
anal=[];
for j=1:n-1
   anal=[anal,sin(2*pi*x)*cos(2*pi*(j/n)^2)];
end

% the maximal difference on the grid 
% between the numerical solution u and the analytic solution anal
FpointNeu1=max(max(abs(anal-u)));


%***************  plotting the solutions **********************
if plot==0
   return
else
   
% order the solution vectors as matrices

mat=[]; % matlab solution u
sol=[]; % analytic solution u0
n=n-1;
for j=0:n-1
   sol=[sol;anal(j*n+1:j*n+n)];
   mat=[mat;u(j*n+1:j*n+n)];
end
n=n+1;

surf(x,x,mat)
shading interp
xlabel('x-direction')
ylabel('y-directon')
title('numerical solution')

figure

surf(x,x,sol)
shading interp
xlabel('x-direction')
ylabel('y-directon')
title('analytic solution')

end %if