function cpr_test ( )

%*****************************************************************************80
%
%% cpr_test() tests cpr().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    20 March 2024
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../cpr' );

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'cpr_test():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  Test cpr(), which uses the Chebyshev Proxy Rootfinder\n' );
  fprintf ( 1, '  to determine all real zeros of a smooth function\n' );
  fprintf ( 1, '  over an interval [a,b].\n' );

  bessel_test ( );
  newton_test ( );
  jenkins_test ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'cpr_test():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  timestamp ( );

  rmpath ( '../cpr' );

  return
end
function bessel_test ( )

%*****************************************************************************80
%
%% bessel_test() tests cpr() on the J0 Bessel function.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    21 March 2024
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'bessel_test():\n' );
  fprintf ( 1, '  Seek real roots of the J0 Bessel function.\n' );
%
%  The search interval:
%
  a = 0.0;
  b = 40.0;
%
%  Chebyshev degree.
%
  N = 50;
%
%  Define the string with the name of the function.
%
  f = @(x) besselj ( 0, x );
%
%  Request roots in interval.
%
  [ froots, Einterstitial, root_residuals ] = cpr ( f, a, b, N );
%
%  Report.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Interstitial interpolation error norm = %g\n', Einterstitial );
  fprintf ( 1, '  Number of roots = %d\n', length ( froots ) );
  fprintf ( 1, '  Roots computed by CPR:\n' );
  disp ( froots );
  fprintf ( 1, '  Maximum residual at roots = %g\n', max ( abs ( root_residuals ) ) );
  exact = [ ...
     2.404825557695773; ...
     5.520078110286311; ...
     8.653727912911012; ...
    11.79153443901428; ...
    14.93091770848779; ...
    18.07106396791092; ...
    21.21163662987926; ...
    24.35247153074930; ...
    27.49347913204025; ...
    30.63460646843198; ...
    33.77582021357357; ...
    36.91709835366404 ]; 
  fprintf ( 1, '  Exact roots:\n' );
  disp ( exact );
  fprintf ( 1, '  Exact f(roots):\n' );
  disp ( f ( exact ) );
  fprintf ( 1, '  f ( froots ):\n' );
  disp ( f ( froots ) );

  return
end
function newton_test ( )

%*****************************************************************************80
%
%% newton_test() tests cpr() on the Newton example.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    21 March 2024
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'newton_test():\n' );
  fprintf ( 1, '  Seek real roots of the Newton example x^3-2x-5=0.\n' );
%
%  The search interval:
%
  a = -5.0;
  b = +5.0;
%
%  Chebyshev degree.
%
  N = 10;
%
%  Define the string with the name of the function.
%
  f = @(x) newton_example ( x );
%
%  Request roots in interval.
%
  [ froots, Einterstitial, root_residuals ] = cpr ( f, a, b, N );
%
%  Report.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Interstitial interpolation error norm = %g\n', Einterstitial );
  fprintf ( 1, '  Number of roots = %d\n', length ( froots ) );
  fprintf ( 1, '  Maximum residual at roots = %g\n', max ( abs ( root_residuals ) ) );
  exact = [ 2.094551481542328 ]; 
  fprintf ( 1, '  Exact real root   : %20.16g\n', exact );
  fprintf ( 1, '  Computed real root: %20.16g\n', froots(1) );
  fprintf ( 1, '  Exact f(roots)    : %20.16g\n', f ( exact ) );
  fprintf ( 1, '  Computed f(roots) : %20.16g\n', f ( froots(1) ) );
  fprintf ( 1, '  Error             : %20.16g\n', abs ( froots(1) - exact ) );

  return
end
function value = newton_example ( x )

%*****************************************************************************80
%
%% newton_example() evaluates the Newton example function.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    21 March 2024
%
%  Author:
%
%    John Burkardt
%
  value = x.^3 - 2.0 * x - 5.0;

  return
end
function jenkins_test ( )

%*****************************************************************************80
%
%% jenkins_test() tests cpr() on the Jenkins polynomial.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    25 March 2024
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'jenkins_test():\n' );
  fprintf ( 1, '  Seek real roots of the Jenkins polynomial\n' );
  fprintf ( 1, '  p(x) = x^4 + 5.6562x^3 + 5.8854x^2\n' );
  fprintf ( 1, '             + 7.3646x   + 6.1354\n' );
%
%  The search interval:
%
  a = -5.0;
  b = 5.0;
%
%  Chebyshev degree.
%
  N = 10;
%
%  Define the string with the name of the function.
%
  f = @(x) jenkins_example ( x );
%
%  Request roots in interval.
%
  [ froots, Einterstitial, root_residuals ] = cpr ( f, a, b, N );
%
%  Report.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Interstitial interpolation error norm = %g\n', Einterstitial );
  fprintf ( 1, '  Maximum residual at roots = %g\n', max ( abs ( root_residuals ) ) );
  fprintf ( 1, '  Number of roots = %d\n', length ( froots ) );
  format ( 'long' );
  fprintf ( 1, '  Computed roots:\n' );
  disp ( froots );
  exact = [ ...
     -4.674054017161709; ...
     -1.0000 ]; 
  fprintf ( 1, '  Exact roots:\n' );
  disp ( exact );
  fprintf ( 1, '  f (Computed roots ):\n' );
  disp ( f ( froots ) );
  fprintf ( 1, '  f(Exact roots):\n' );
  disp ( f ( exact ) );


  return
end
function value = jenkins_example ( x )

%*****************************************************************************80
%
%% jenkins_example() evaluates the Jenkins example function.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    25 March 2024
%
%  Author:
%
%    John Burkardt
%
  value = x.^4 + 5.6562*x.^3 + 5.8854*x.^2+7.3646*x+6.1354;

  return
end
function timestamp ( )

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

  return
end
