function NpointDi1=NpointDi1(n,plot)

% calculating the numerical solution of

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

% 9point-star
% Dirichlet 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 Dirichlet boundary condition g ***************

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

% lower boundary go: y=0
go=sin(2*pi*x);

% left boundary  gl: x=0
gl=0*x;

% right boundary gr: x=1
gr=gl;

% upper boundary gu: y=1
gu=go;

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

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

% to get convergence of the order 4, it is necessary to apply a 5point star on the right side p
% this will be done be applying matrix B on the right side

% B consists of submatrices L
e=ones(n-1,1);
L=spdiags([ e e 8*e e e],[ 1 n+1 n+2 n+3 2*n+3],n-1,3*(n+1));
L=1/12*L;

B=L;
for j=1:(n-2)
   B=[B,spalloc(j*(n-1),n+1,0)
      spalloc(n-1,j*(n+1),0),L];
end

% p will be the right side of the finite difference equation A*u=p
p=(B*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^2*1/6*(go(1:n-1)   +4*go(2:n)  +go(3:n+1));

% left border
p(1:n-1:(n-1)*(n-2)+1)  =p(1:n-1:(n-1)*(n-2)+1)  +n^2*1/6*(gl(1:n-1)   +4*gl(2:n)  +gl(3:n+1));

% right border
p(n-1:n-1:(n-1)^2)      =p(n-1:n-1:(n-1)^2)      +n^2*1/6*(gr(1:n-1)   +4*gr(2:n)  +gr(3:n+1));

% upper border
p((n-1)*(n-2)+1:(n-1)^2)=p((n-1)*(n-2)+1:(n-1)^2)+n^2*1/6*(gu(1:n-1)   +4*gu(2:n)  +gu(3:n+1));

% each corner is member of to borders, that's why it is necessary to undo one addition

% down left
p(1)            =p(1)            -n^2/6*go(1);

% down right
p(n-1)          =p(n-1)          -n^2/6*gr(1);

% upper left
p((n-1)*(n-2)+1)=p((n-1)*(n-2)+1)-n^2/6*gl(n+1);

% upper right
p((n-1)^2)      =p((n-1)^2)       -n^2/6*gu(n+1);

p=p';
%***************  Matrix A **********************

e=ones((n-1)^2,1);
A=spdiags([-e -4*e -e -4*e 20*e -4*e -e -4*e -e],[-n,-n+1,-n+2,-1,0,1,n-2,n-1,n],(n-1)^2,(n-1)^2);

% every time, when an entry of the star touches the border, it is necessary to put its value to zero

for j=1:n-2
   A(j*(n-1),j*(n-1)+1)=0; %right   1 n-2
   A(j*(n-1)+1,j*(n-1))=0; %left    1 n-2
end
for j=1:n-1
   A(j*(n-1),(j-1)*(n-1)+1)=0; %up right   1 n-1
   A((j-1)*(n-1)+1,j*(n-1))=0; %down left  1 n-1
end
for j=1:n-3
   A(j*(n-1),(j+1)*(n-1)+1)=0; %down right  1 n-3
   A((j+1)*(n-1)+1,j*(n-1))=0; %up left     1 n-3 
end
A=n^2*1/6*A;

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

u=A\p;
u=u';


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

% the analytic solution as a vector

x=linspace(1/n,1-1/n,n-1);
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
NpointDi1=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([0:1/n:1],[0:1/n:1],[go;gl(2:n)',mat,gr(2:n)';gu])
shading interp
xlabel('x-direction')
ylabel('y-directon')
title('numerical solution')

figure

surf([0:1/n:1],[0:1/n:1],[go;gl(2:n)',sol,gr(2:n)';gu])
shading interp
xlabel('x-direction')
ylabel('y-directon')
title('analytic solution')

end %if