%% BVP_JACOBI solves a two point boundary value problem using the Jacobi method.
%
%  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 );
  a(1,1) = 1.0;
  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
  a(n,n) = 1.0;
%
%  The right hand side.
%
  b = zeros ( n, 1 );
  b(1) = ua;
  b(2:n-1) = 0.0;
  b(n) = ub;
%
%  Jacobi iteration, from x = 0.
%
  u = zeros ( n, 1 );
  unew = zeros ( n, 1 );

  figure ( 1 )
  clf
  hold on

  for it = 1 : 100
    unew(1) = b(1);
    unew(2:n-1) = ( dx^2 * b(2:n-1) + u(1:n-2) + u(3:n) ) / 2.0;
    unew(n) = b(n);
    r = norm ( a * u - b );
    fprintf ( 1, '  %d  %g\n', it, r );
    u(1:n) = unew(1:n)

    plot ( x, u, 'b-', 'Linewidth', 3  )
    pause

  end

  hold off


