function h = plot_Fneq_panel(constraints, universal, timesteps, waters)

max_timestep = 7.0;
ntimesteps = length(timesteps);
nwaters = length(waters);
colors = jet(ntimesteps);

dt_ni = zeros(nwaters, ntimesteps);
Fneq_ni = zeros(nwaters, ntimesteps);
dFneq_ni = zeros(nwaters, ntimesteps);

hold on;
for i = 1:ntimesteps
  try
    for j = 1:nwaters
      timestep = timesteps{i};
      nwat = waters{j};
      filename = sprintf('analysis3/TIP3P-%d-%s-timestep-%04.2f.out', nwat, constraints, timestep)
      x = textread(filename);
      dt_ni(j,i) = timestep;
      
      %index = 9; % which ntimesteps to use (256)
      %index = 10; % which ntimesteps to use (512)
      index = 11; % which ntimesteps to use (1024)
      %index = size(x,1); % which ntimesteps to use
      Fneq_ni(j,i) = x(index,9);
      dFneq_ni(j,i) = x(index,10);

      if universal
	ndof = nwat*9;
	if strcmp(constraints, 'constrained')
	  ndof = nwat*(9-3);
	end
	Fneq_ni(j,i) = Fneq_ni(j,i) / nwat;
	dFneq_ni(j,i) = dFneq_ni(j,i) / nwat;
      end
    end
  catch
    % Could not read file---the timestep stabiity limit was reached.
    disp(sprintf('stability limit reached at i = %d', i));
    disp(filename);
    ntimesteps = i-1
    dt_ni = dt_ni(:, 1:ntimesteps);
    Fneq_ni = Fneq_ni(:, 1:ntimesteps);
    dFneq_ni = dFneq_ni(:, 1:ntimesteps);
    max_timestep = timesteps{ntimesteps};
    break
  end
end

h = errorbar(dt_ni', Fneq_ni', 2*dFneq_ni', '.-');
set(h, 'markersize', 10);

return
