function graphHarmonic(nValues, trials)
  % Graphs the theoretical harmonic numbers versus the observed peaks.
  % USAGE:   graphHarmonic(nValues, trials)
  %          nValues is a vector of positive values for n
  %          trials is number of trials for each n Value
  %             (trials is 1000 by default)
  if nargin < 2
      trials = 1000;
  end
  largest = max(nValues);
  fractions = 1 ./ (1:largest);
  theoretical = cumsum(fractions);
  hold off;
  plot(theoretical,'r');
  xlabel('value of n');
  ylabel('number of peaks');
  hold on;
  
  % Run trials using the harmonicTest function
  for n = nValues
      avgPeaks = harmonicTest(n, trials);
      plot(n,avgPeaks,'bo');
  end
  
  hold off;