  m = 5.0;
  bVals = 0.1 : 0.05 : 5;
  kVals = 1.5 : 0.05 : 5;

  [ kGrid, bGrid ] = meshgrid ( bVals, kVals );

  peakVals = nan ( size ( kGrid ) );

  tic;

  parfor ij = 1 : numel(kGrid)

    [ T, Y ] = ode45 ( @(t,y) ode_system ( t, y, m, bGrid(ij), kGrid(ij) ), ...
      [0, 25],  [0, 1] );

    peakVals(ij) = max ( Y(:,1) );

  end

  toc;
