
function [u,res]=FEM1DPureNeu17(f,a,b,g1,g2,x=linspace(a,b,40+1))
%FEM 1D solver for -u''=f on (a,b) 
%with the puere Neumann condition 
%on bnd: u'(a)=g1, u'(b)=g2 
%INPUT
%f - function handle for f rhs function (should work for vectors)
%a,b - interval endswitch
%g1,g2 - Neumann bnd alues at a and b
%x- vector with mesh points; a column!!!!
%OUTPUT
%u the vector of values at nodal points s.t. u(1)=0 i.e. equal t ozero at a
%res=norm(A*u-G,'inf') - max norm of residual 
%
%EXAMPLE
%N=50;a=0;b=2;x=linspace(a,b,N+1)';f=@(x) ones(size(x));g1=0;g2=-2;
%[u,res]=FEM1DPureNeu17(f,a,b,g1,g2,x);norm(u+0.5.*x.*x,'inf')
%we know solution here : u=-0.5*x.*x
%different mesh
%x=0.5*(a+b)-0.5*(b-a)*cos(pi*(linspace(a,b,N+1)-a)/(b-a));x=x';
%[u,res]=FEM1DPureNeu17(f,a,b,g1,g2,x);norm(u+0.5.*x.*x,'inf')

A=FEM1Dmats(x');

N=length(x)-1; 
hh=zeros(1,N); 
#hh(k)=x(k+1)-x(k)
hh=x(2:N+1)-x(1:N);
h=max(hh);

F=f(x);
G=zeros(size(F));
G(1:N)=0.5*F(1:N).*hh;
G(2:N+1)+=0.5*F(2:N+1).*hh;
G(1)+=-g1;
G(N+1)+=g2;

if(abs(sum(G))>1e-12)
 printf("f and g1 and g2 do NOT satisfy the compability solution\n");
 printf("Terminating the solver\n");
 return; 
end 

A1=A(2:N+1,2:N+1);
G1=G(2:N+1);
u=A1\G1;
u=[0;u];
%we compute the max norm of the residual G-A*u
res=norm(G-A*u,'inf');
end






