function sol=femsolver(x,lbcparam,rbcparam,afu,ffu,xffu)
%% Does the dirty work for the master script.
%% The arguments are explained there.
lbctype=lbcparam(1); a0=lbcparam(2); b0=lbcparam(3); u0=lbcparam(4);
rbctype=rbcparam(1); a1=rbcparam(2); b1=rbcparam(3); u1=rbcparam(4);

%% Generate the lengths of intervals between x-nodes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m=length(x)-2;

for i=1:m+1,
  h(i)=x(i+1)-x(i);
end

h=h';

%% Forming the rhs "load vector"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Middle components as usual
for i=1:m;
  b(i+1)=(quad8(xffu,x(i),x(i+1)) ...
              -x(i)*quad8(ffu,x(i),x(i+1)))/h(i);
  b(i+1)=b(i+1)+(x(i+2)*quad8(ffu,x(i+1),x(i+2)) ...
              - quad8(xffu,x(i+1),x(i+2)))/h(i+1);
end;

%% First component depending on lbctype
if lbctype == 1
  b(1)=(x(2)*quad8(ffu,x(1),x(2)) ...
  -quad8(xffu,x(1),x(2)))/h(1)+b0*aa(0);
else
  b(1)=u0;
end

%% Last component depending on rbctype
if rbctype == 1
  b(m+2)=(quad8(xffu,x(m),x(m+1)) ... 
  -x(m)*quad8(ffu,x(m),x(m+1)))/h(m)+b1*aa(1);
else
  b(m+2)=u1;
end 

b=b';

%% Forming the lhs "stiffness matrix"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
%% Main diagonal
for i=1:m,
  amat(i+1,i+1)=quad8(afu,x(i),x(i+1))/h(i)^2 ...
                +quad8(afu,x(i+1),x(i+2))/h(i+1)^2;
end

%% Side diagonals
for i=1:m+1,
  amat(i,i+1)=-quad8(afu,x(i),x(i+1))/h(i)^2;
  amat(i+1,i)=amat(i,i+1);
end

%% First row depending on lbctype
if lbctype == 1
  amat(1,1)=a0*aa(0)+quad8(afu,x(1),x(2))/h(1)^2;
else
  amat(1,1)=1;   
  amat(1,2)=0;
end

%% Last row depending on rbctype
if rbctype == 1
  amat(m+2,m+2)=a1*aa(1)+quad8(afu,x(m+1),x(m+2))/h(m+1)^2;
else
  amat(m+2,m+2)=1;
  amat(m+2,m+1)=0;
end

%% Solution of the resulting system (amat*sol=b)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sol=amat\b;


