


function [u,hh,A,B]=FEM1Dsolver(f,x,al=0,be=0,c=0)
%Solver of -u''+cu=f u(a)=al u(b)=be  
%INPUT
%f function handle to the rhs function 
%x the vector with the mesh a=x(1) b=x(end)
%al,be - the Dirichelt bnd  values i.e. u(a)=al u(b)=be
%c- the reaction coefficient (default =0)
%OUTPUT
%u the solution vector ui.e. u(k) approximates solution at x(k)
%hh the mesh vector (hh(k)=x(k+1)-x(k))
%A,B - the stiffness and the mass matrices of linear FEM on the mesh x  
%EXAMPLE:
%a=0;b=pi;c=1;al=be=0;f=@(x) (1+c).*sin(x);x=linspace(a,b,20)';[u]=FEM1Dsolver(f,x,al,be,c);
%a=0;b=1;c=1;al=sin(a);be=sin(b);f=@(x) (1+c).*sin(x);x=linspace(0,1,20)';x=a+(b-a)*sin(x)/sin(1);
%[u]=FEM1Dsolver(f,x,al,be,c);


  [A,B,h,hh]=FEM1Dmats(x);
  N=length(hh);
  %Since we have dirichlet bnd cond we cat extract submatrices
  iN=2:N;
  iA=A(iN,iN);iB=B(iN,iN);
  ix=x(iN);
  %FEM 1D RHS 
  F=f(ix).*(hh(1:N-1)+hh(2:N))*0.5;
  F(1)+=-A(2,1)*al;
  F(N-1)+=-A(N,N+1)*be;
  u=(iA+c*iB)\(F); %solves (iA+c*iB)*u=F
  u=[al;u;be];
end




function [A,B,h,hh]=FEM1Dmats(x,N=20)
%FEM creates matrises stiff and mass for  -u''+c*u=f 
%INPUT
%x either a vector od nodal points a=x(1)<...<x(N+1)=b
%  or a vector of length 2 [a,b] - the ends of the interval
%N (default =20) the no of equidistant nodal points minus one - in case of 
%  x being the 2 element vector with the ends of the interval 
%  (otherwise this parameter is ignored) 
%OUTPUT
%A -stiffness matrix A=(\int_a^b \phi_k'\phi_l'dx)_{k,l=0}^N
%B mass matrix B=(\int_a^b \phi_k\phi_ldx)_{k,l=0}^N
%h -=max(hh) maximal size of subinterval
%hh - vector of the sizes of subintervals i.e. hh(k)=x(k+1)-x(k)
%EX: 
%a=0;b=1; [A,B,h]=FEM1Dmats([a,b],10);h*full(A)
%a=0;b=1; x=linspace(0,1,11)';x=a+(b-a)*sin(x)/sin(1); [A,B]=FEM1Dmats(x);

if(rows(x)==1) x=x';
end
if(length(x)==2)
 a=x(1);b=x(2);
 x=linspace(a,b,N+1)';
endif

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);

mdA=1./hh(1:N-1)+1./hh(2:N);
mdA=[1/hh(1);mdA;1/hh(N)]; 
mdB=hh(1:N-1)+hh(2:N); 
mdB=[hh(1);mdB;hh(N)];
subd=supd=hh;
#stiffness
A=sparse(diag(-1.0./subd,-1)+ diag(mdA) + diag(-1.0./supd,1));
#mass
B=sparse(diag(subd./6,-1)+ diag(mdB/3) + diag(supd/6,1));

end
