function peakVals = spring_sweep ( bVals, kVals )

%*****************************************************************************80
%
%% spring_sweep() finds peak values of spring_ode() over a mesh of parameter values.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    21 March 2011
%
%  Author:
%
%    John Burkardt.
%
%  Input:
%
%    real bVals(M), a vector of B values.
%
%    real kVals(N), a vector of K values.
%
%  Output:
%
%    real peakVals(M,N), a vector containing the maximum value of
%    the ODE solution over a fixed time interval, for every combination
%    of B and K values.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'spring_sweep():\n' );
  fprintf ( 1, '  Tabulate the maximum solution value of the spring\n' );
  fprintf ( 1, '  ODE over a range of K and B values.\n' );
%
%  Form a grid of all pairs of K and B:
%
  [ kGrid, bGrid ] = meshgrid ( bVals, kVals );
%
%  Define an array to hold the results, and initialize it to NAN.
%
  peakVals = nan ( size ( kGrid ) );
%
%  Solve the ODE for every pair of K and B values.
%
  [ m, b, k, t0, y0, tstop ] = spring_parameters ( );

  tspan = [ t0, tstop ];
  f = @spring_deriv;

  for ij = 1 : numel(kGrid)
%
%  Replace default vallues of B and K.
%
    b = bGrid(ij);
    k = kGrid(ij);
    [ m, b, k, t0, y0, tstop ] = spring_parameters ( m, b, k, t0, y0, tstop );
%
%  Solve the ODE.
%
    [ t, y ] = ode45 ( f, tspan, y0 );
%
%  Retrieve the maximum value achieved by this solution.
%
    peakVals(ij) = max ( y(:,1) );

  end

  return
end
