

function [A,B,h,hh]=FEM1Dmats(x,N=20)
%FEM creates matrises stiff and mass for  -u''+c*u=f 
%
%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=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
