ellfem.m

various matlab programs to slove various problems

function [uh,x] = ellfem(Nx,I,alpha,beta,gamma,f,bc,stabfun)
a = I(1); 
b = I(2); 
h = (b-a)/Nx; 
x = [a+h/2:h:b-h/2];
alpha = eval(alpha); 
beta = eval(beta); 
gamma = eval(gamma);
f     = eval(f);     
rhs  = 0.5*h*(f(1:Nx-1)+f(2:Nx));    
if nargin == 8
  [Afe,rhsbc] = femmatr(Nx,h,alpha,beta,gamma,stabfun);
else
  [Afe,rhsbc] = femmatr(Nx,h,alpha,beta,gamma);
end
[L,U,P]   = lu(Afe);
rhs(1)    = rhs(1)-bc(1)*(-alpha(1)/h-beta(1)/2+h*gamma(1)/3+rhsbc(1));
rhs(Nx-1) = rhs(Nx-1)-bc(2)*(-alpha(Nx)/h+beta(Nx)/2+h*gamma(Nx)/3+rhsbc(2));
rhs = P*rhs'; 
z = L \ rhs; 
w = U \ z; 
uh = [bc(1), w', bc(2)]; 
x  = [a:h:b];
return