%% BVP_GAUSS solves a two point boundary value problem with Gauss elimination.
%
%  The physical region [AX,BX]:
%
  ax = 5.0;
  bx = 10.0;
%
%  The boundary conditions:
%
  ua = 100.0;
  ub = 50.0;
  uap = -10.0;
  ubp = -10.0;
%
%  Set up the nodes.
%
  n = 21;
  x = linspace ( ax, bx, n );
  dx = ( bx - ax ) / ( n - 1 );
%
%  Set up the matrix:
%  We are approximating -uxx, so the signs in the finite difference
%  approximation of d2udx2 must be switched here.
%
  a = zeros ( n, n );
  for i = 2 : n - 1
    a(i,i-1) = -1.0 / dx^2;
    a(i,i  ) =  2.0 / dx^2;
    a(i,i+1) = -1.0 / dx^2;
  end
%
%  The right hand side.
%
  b = zeros ( n, 1 );
  b(2:n-1) = 0.0;
%
%  The solution of A*x=b is 0 if b is zero...
%  Oh, we forgot the boundary conditions!
%
  a(1,1) = 1.0;
  b(1) = ua;
  a(n,n) = 1.0;
  b(n) = ub;

  u = a \ b;

  clf
  hold on 
  plot ( x, u, 'b-', 'Linewidth', 3  )
  grid on

  pause ( )

  for i = 2 : n - 1
%
%  Fix the sign.
%
%   b(i) = 1.0 / ( 1.0 + ( x(i) - 9.0 )^2 );
%
%  Make bump at 9 bigger.
%
%   b(i) = 1.0 / ( 1.0E-01 + ( x(i) - 9.0 )^2 );
%
%  More bigger.
%
%   b(i) = 1.0 / ( 1.0E-02 + ( x(i) - 9.0 )^2 );
%
%  Most bigger.
%
%   b(i) = 1.0 / ( 1.0E-03 + ( x(i) - 9.0 )^2 );
%
%  Spread bump out to see response is not a piecewise linear.
%
%   b(i) = 10.0 / ( 1.0E-01 + ( x(i) - 9.0 )^2 );
%
%  Heat and cool.
%
    b(i) = 1 * x(i).^2 * sin ( pi * x(i) );

  end

  u = a \ b;

  plot ( x, u, 'ro-', 'Linewidth', 3 )
%
%  Plot the source term
%
  plot ( x(2:n-1), b(2:n-1), 'g-', 'Linewidth', 3 )

  hold off

  pause
%
%  Solve problem with left boundary condition u'(a) = +10.
%
  a(1,1) = - 1.0 / dx;
  a(1,2) = + 1.0 / dx;
  b(1) = + 10.0;
  u = a \ b;

  clf
  hold on
  plot ( x, u, 'bo-', 'Linewidth', 3 )
  plot ( x(2:n-1), b(2:n-1), 'g-', 'Linewidth', 3 )
  grid on
  hold off
