function fem_neumann_test ( )

%*****************************************************************************80
%
%% fem_neumann_test() tests fem_neumann()..
%
%  Discussion:
%
%    This script sets up and runs a finite element simulation code 
%    for a 1D reaction/diffusion equation.
%
%    The dynamics are given by the following reaction/diffusion equation
%    for the function W(T,X):
%
%      W_t = W_xx + NL(W,c),
%
%    where NL(W,c) is a polynomial in W:
%
%       NL(W,c) = c(1) + c(2) * W + c(3) * W^2 + c(4) * W^3
%
%    with zero Neumann boundary conditions at X = 0.0 and X = 1.0:
%
%      W_x(T,0.0) = 0.0
%      W_x(T,1.0) = 0.0
%
%    and initial condition at T = 0.0:
%
%      W(0,X) = sin ( pi * X ).
%
%    The problem is to be solved for 0.0 <= T <= 4.0, 0.0 <= X <= 1.0.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    22 January 2019
%
%  Author:
%
%   Eugene Cliff
%
%  Reference:
%
%    Jeffrey Borggaard, John Burkardt, John Burns, Eugene Cliff,
%    Working Notes on a Reaction Diffusion Model: a Finite Element Formulation.
%
  addpath ( '../fem_neumann' );

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'fem_neumann_test():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  fem_neumann() solves a 1D time-dependent reaction/diffusion\n' );
  fprintf ( 1, '  problem with a nonlinear term and Neumann boundary conditions.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The equation is discretized using the finite element method.\n' );
%
%  Set initial condition.
%  The first option uses a centered hat.
%
  if ( false )
    w_0 = @(x) 1.0 + basic_hat ( 8.0 * x - 4.0 );
  else
    w_0 = @(x) sin(pi*x);
  end
%
%  Set the polynomial coefficients of the reaction term.
%  Here, c = 0*x -1*x + 0*x^2 + 1*x^3 = x * ( x^2 - 1 )
%
  c = [ 0.0; -1.0; 0.0; 1.0 ];
%
%  Select 101 times for output.
%
  nt = 101;
  t = linspace ( 0.0, 4.0, nt );
%
%  The number of spatial grid points is N+1;
%
  n = 32;
  x = linspace ( 0.0, 1.0, n + 1 );
%
%  T is a vector of 101 times.
%  W is an array of 101x33 function values.
%
  [ T, W ] = rd_lin_spline ( w_0, t, n, c );
%
%  Plot the results.
%
  plot_rd ( t, x, W, c )
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'fem_neumann_test:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../fem_neumann' );

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% timestamp() prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

