function pvals = ode_sweep_parfor ( bvals, kvals )

%*****************************************************************************80
%
%% ode_sweep_parfor() finds peak values of an ODE over a mesh of parameter values.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    18 December 2023
%
%  Author:
%
%    John Burkardt.
%
%  Input:
%
%    real bvals(M), a vector of B values.
%
%    real kvals(N), a vector of K values.
%
%  Output:
%
%    real pvals(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.
%

%
%  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.
%
  pvals = nan ( size ( kGrid ) );
%
%  Solve the ODE for every pair of K and B values.  (M is fixed at 5.)
%
  m = 5.0;

  parfor ij = 1 : numel(kGrid)
%
%  Solve the ODE over the time interval 0 <= T <= 25, with
%  initial conditions X(0) = 0, X'(0) = 1.
%
    [ T, Y ] = ode45 ( @(t,y) ode_sweep_system ( t, y, m, bGrid(ij), kGrid(ij) ), ...
      [0, 25],  [0, 1] );
%
%  Retrieve the maximum value achieved by this solution.
%
    pvals(ij) = max ( Y(:,1) );

  end

  return
end
