
function [u,gr,h,A,B,C,S,F]=FEM1Dsolver(f,a,b,x=[],al=0,be=0,c=0,ad=0,di=1.0)
#FEM solver for -di*u''+ad*u'+ c*u=f with
#Neumann/Robin bnd cond -u'(a)+cau(a)=al u'(b)+cb*u(b)=be
#or Dirichlet bnd : u(a)=al u(b)=be
#Input: f - function handle
#a,b- ends of interval
#x- no of mesh points OR a vector with mesh points : e.g. (xk=a+k*h) - default such with h=(b-a)/20 if x=[] then mesh equdistant
#al, beta - rhs of  bnd conditions
#c - either a value of c (Dirichelt bnd) or a vector [c,ca,cb] - then we have Robin bnd cond
#ad - advection parameter
#di - diffusion parameter 
#output:
#u - FEM1D linear  solution on these mesh - the rhs computed by trapezoidal rule
#gr =- grid
#h - max mesh parameter i.e. h=max h_k h_k=x_{k}-x_{k-1} k=1,..,N
#A,B,C,S - matrices - diffusion, mass, advection - final FEM matrix
#F - rhs of the FEM system 
#Ex:
#[u,t]=FEM1Dsolver(@sin,0,pi);plot(t,u)
#[u,t]=FEM1Dsolver(@sin,a=-1,b=2,[],sin(a),sin(b));plot(t,u) #nonhomogenous D bc
#[u,t]=FEM1Dsolver(@sin,a=-1,b=2,[],0.5*sin(a),0.5*sin(b),1.0);plot(t,u,t,sin(t)/2)
#[u,t]=FEM1Dsolver(@(x)2*sin(x)+cos(x),a=-1,b=2,[],sin(a),sin(b),1.0,1.0);plot(t,u,t,sin(t))
#[u,t]=FEM1Dsolver(@sin,0,pi,[],-cos(0),cos(pi),[0,0,0]);plot(t,u)#Neumanna bnd
#[u,t]=FEM1Dsolver(@sin,a=-2,b=1,[],-cos(a)+sin(a),cos(b)+sin(b),[0,1,1]);plot(t,u)#Robin bnd
#[u,t]=FEM1Dsolver(@(x) 1,0,1,100,0,-1,[0,0,0]);plot(t,u)#neumann  bnd 
#[u,t]=FEM1Dsolver(@(x) 1,0,1,100,0,0,[0,0,0]);# neumann bnd cond -  WRONG BND cond
#[u,t]=FEM1Dsolver(@(x) 100.0,0,2,100,0,1);plot(t,u) #very concave -u''=100
#[u,t]=FEM1Dsolver(@(x) -100.0,0,2,100,0,1);plot(t,u) #very convex -u''=-100
#[u,t]=FEM1Dsolver(@(x) (100+10*x)*(0.5-(x>0)),-1,2,100,-1,1);plot(t,u) #pu
#[u,t,h,A,B,S,F]=FEM1Dsolver(@sin,0,pi);plot(t,u)
#denser grid on ((b-a)/2,b] subinterval
#x1=linspace(0,pi/2,20);x2=linspace(pi/2,pi,1000);x=[x1,x2(2:999)];[u,t]=FEM1Dsolver(@sin,0,pi,x);plot(t,u-sin(t),";err;");
#Advection 
#eps=1e-1;f=@(x) (abs(x)>0.1)./eps;[u,t]=FEM1Dsolver(f,-1,1,100,0,0,0,ad=0);plot(t,u)
#eps=1e-1;f=@(x) (abs(x)>0.1)./eps;[u,t]=FEM1Dsolver(f,-1,1,100,0,0,0,ad=10);plot(t,u)
#eps=1e-1;f=@(x) (abs(x)>0.1)./eps;[u,t]=FEM1Dsolver(f,-1,1,100,0,0,0,ad=-10);plot(t,u)
#very small diffusion with respect to advection
#eps=1e-1;f=@(x) d*(abs(x)>0)./(eps);[u,t]=FEM1Dsolver(f,-1,1,100,0,0,0,-1,1e-3);plot(t,u)#ti does not work. Too small h... with respect to ad/di


#Code starts
type=0; #0 Dirichlet; 1 Robin; 2 neumann


if(isempty(x))
 N=round((b-a)*20);
 x=linspace(a,b,N+1);
else 
   if(length(x)==1)
     N=x; 
     x=linspace(a,b,1+N);
   else
     N=length(x)-1;
     a=x(1);
     b=x(N+1);
   endif  
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));
#\int u'v advection local (\phi_k,\phi_{k+1})=0.5 [-1,]
supd=ones(1,N);subd=-supd;
C=sparse(diag(0.5*subd,-1)+diag(0.5*supd,1));
C(1,1)=-0.5;C(N+1,N+1)=0.5;

F=zeros(N+1,1);
for k=0:N,
 F(k+1)=f(x(k+1));
endfor
#scaling
F=0.5*F.*mdB';
S=di*A+ad*C+c(1)*B;

if(length(c)==1)
 #Dirichlet bnd 
 F(2)-=S(2,1)*al; #-al*\int_a^{x_1} \phi_0'\phi_1'=al/h_0
 F(N)-=S(N,N+1)*be;
 S=S(2:N,2:N);
 F=F(2:N);
else 
 type=1;
#Robin bnd weak formul: 
#\int u'v' + cuvdt + u'(a)v(a)-u'(b)v(b)=\inf fv 
#\int..dt + cau(a)v(a)+cbu(b)v(b)=\intfv +alv(a)+ be v(b)
 F(1)+=al;F(N+1)+=be; #adding contribution of  alv(a)+ be v(b)
 S(1,1)+=c(2);
 S(N+1,N+1)+=c(3); # ca*u(a)v(a)+cb*u(b)v(b) to stiffness matrix
 if(norm([ad,c(2:3)],'inf')<1e-19) 
   type=2;
  ind=1:N;
   #case of pure Neumann condition we assume that u(N+1)=u(b)=0 to get a unique solution
   S=S(ind,ind);
   if(sum(F)>h)
     printf("Pure Neumann cond sum(F(k)<>0! Problem has no solution!\n");
   endif
   F=F(ind);
 endif
endif

u=S\F;
if(type==0)
 u=[al; u; be];
 else 
 if(type==2)
  #u coefficients of the solution in the basis (\phi_0,..\phi_{N-1}, 1)
  #here 
  u=[u;0];
 endif
endif
gr=x'; #returns grid 

endfunction

