function [xv,uv,stime] = fdbvp_2nd(n,a,b,ua,ub,sigma,f)
%
h = (b-a)/n;
xv = a:h:b; xv=xv';
fv = zeros(n+1,1);
sigmav = zeros(n+1,1);

for i = 1:n+1
    x = xv(i);
    sigmav(i) = eval(sigma);
    fv(i) = eval(f);
end



av = (2*ones(n-1,1)+h*h*sigmav(2:n))/(h*h);
bv = -ones(n-1,1)/(h*h); bv(1)=0;
cv = -ones(n-1,1)/(h*h); cv(n-1)=0;

fv = fv(2:n); 
fv(1) = fv(1) + ua/(h*h); fv(n-1) = fv(n-1) + ub/(h*h);

% general solver  
%
% A = (1/(h*h))*sparse(diag(2*ones(n-1,1)+h*h*sigmav(2:n),0) ...
%                    - diag(ones(n-2,1),-1) - diag(ones(n-2,1),1));
% uv = A\fv;  
%

tic;
[alpha,beta] = trilu(av,bv,cv);
uv = trilu_solve(alpha,beta,cv,fv);
stime = toc;

uv = [ua;uv;ub];


