
function [u,hh,A,B]=FEM1dDirSol16(x,f,c=0,al=0,be=0)
%solves -u''+cu=f on [a,b] u(a)=0=u(b) by 1d linear  fem;  
%in case x=[a,b] mesh is equidistant on N equal subintervals
%INPUT
%x -mesh in increasing order
%f -function handle for the rhs
%c- nonnegative coeff
%OUTPUT
%u solution in nodal basis 
%h - max length of subintervals
%we create matrices for the full basis (\phi_k)_{k=0}^N
%EXAMPLE
%N=20;f=s=@sin;x=linspace(0,pi,N)';[u,h,A,B]=FEM1dDirSol16(x,f);er=u-s(x);eH1=sqrt(er'*A*er);eL2=sqrt(er'*B*er);[eH1,eL2]
%N=10;f=s=@sin;x1=linspace(0,b=pi/2,N)';x2=linspace(b,pi,2*N)';x=[x1;x2(2:end)];[u,h,A,B]=FEM1dDirSol16(x,f);er=u-s(x);eH1=sqrt(er'*A*er);eL2=sqrt(er'*B*er);[eH1,eL2]

if(columns(x) < rows(x))
 x=x'; %x must be horizontal
end
[A,B,h]=FEM1Dmats(x);
N=length(x)-1;
hh=x(2:N+1)-x(1:N);hh=hh';
F=f(x)'; 
AA=A(2:N,2:N);
BB=B(2:N,2:N);
F=F(2:N);
F=0.5*F.*(hh(1:N-1)+hh(2:N));%scaling F(i)=f(x(i))*(hh(i-1)+hh(i))
ub=zeros(N+1,1);
ub(1)=al;ub(N+1)=be;
Aub=A*ub+c*B*ub;
ui=(AA+c*BB)\(F-Aub(2:N));%solving
u=[al;ui;be]; %adding zero Dirichlet bnd 
end




