%   Author: Reija Ristilä <reija@hut.fi>
%   Last modified: 1999-04-19 14:55:37 EEST

function [u,A,b]=h53b(x,bc)

% Problem: -(a(x)u'(x))' = f(x), 0<x<1
% Files: 
%  - a(x) in file ax.m
%  - f(x) in file fx.m
%  - x*f(x) in file xfx.m
% Arguments: x & bc
%  - for x include beginning and end points (0 and 1); column vector
%  - boundaries: first indicate type: 0=Dirichlet, 1=Neumann,
%                then give appropriate values for both ends:
%                for example bc=[ 0 1 u0 g1];

% function y=ax(x)
% y=1;
% function y=fx(x)
% y=-exp(x-1);
% function y=xfx(x)
% y=x.*fx(x);

if size(x,1) < size(x,2)
  x=x';     % column vector
end
m = size(x,1);
  
%%%% mesh
h=diff(x);

%%%% stiffness matrix A
A=sparse(m,m);
for i=2:(m-1)
  A(i,i)=quad8('ax',x(i-1),x(i))/h(i-1)^2+quad8('ax',x(i),x(i+1))/h(i)^2;
  A(i,i+1)=-quad8('ax',x(i),x(i+1))/h(i)^2;
  A(i+1,i)=A(i,i+1);
end

%%%% load vector b
b=zeros(m,1);
% inner nodes, i=2:m-1
for i=2:(m-1);
  b(i)=(quad8('xfx',x(i-1),x(i))-x(i-1)*quad8('fx',x(i-1),x(i)))/h(i-1);
  b(i)=b(i)+(x(i+1)*quad8('fx',x(i),x(i+1))-quad8('xfx',x(i),x(i+1)))/h(i);
end;

%%%% boundary at x=0, i=1:
if bc(1)==1       % Neumann
  A(1,1)= quad8('ax',x(1),x(2))/h(1)^2;
  b(1)=(x(2)*quad8('fx',x(1),x(2))-quad8('xfx',x(1),x(2)))/h(1)-bc(3)*ax(0);
else              % Dirichlet
  A(1,1)=1;
  b(1)=bc(3);
end

%%%% boundary at x=1, i=m:
if bc(2)==1       % Neumann
  A(m,m)= quad8('ax',x(m-1),x(m))/h(m-1)^2;
  b(m)=(quad8('xfx',x(m-1),x(m))-x(m-1)*quad8('fx',x(m-1),x(m)))/h(m-1)+bc(4)*ax(m);
else              % Dirichlet
  A(m,m)=1; A(m,m-1)=0;
  b(m)=bc(4);
end

u=A\b;
