x3min =  0.5;
x3max =  4;
y3min =  0;
y3max =  3.5;

numSteps = 201;

aliasFact = 10;

stopR = 5;

initvx1 = 0;
initvx2 = 0;
initvx3 = 0.1;

initvy1 = 0;
initvy2 = 0;
initvy3 = 0;

initx1 = 0;
initx2 = 1;
inity1 = 0;
inity2 = 0;

x3arr = linspace(x3min, x3max, numSteps);
y3arr = linspace(y3min, y3max, numSteps);

finalLoner = 0;
for i = 1 : numSteps
  initx3 = x3arr(i);
  for j = 1 : numSteps
    inity3 = y3arr(j);
    try
      sim('ThreeBody', 500);
      
      r12 = sqrt((x1(end) - x2(end)) .^ 2 + (y1(end) - y2(end)) .^ 2);
      r13 = sqrt((x1(end) - x3(end)) .^ 2 + (y1(end) - y3(end)) .^ 2);
      r32 = sqrt((x3(end) - x2(end)) .^ 2 + (y3(end) - y2(end)) .^ 2);
      
      if ((r12 < r13) & (r12 < r32))
        finalLoner(i,j) = 3;
      elseif ((r13 < r12) & (r13 < r32))
        finalLoner(i,j) = 2;
      else
        finalLoner(i,j) = 1;
      end;
      catch,
        finalLoner(i,j) = NaN;
      end;
      
    disp((i - 1) / numSteps + j / numSteps ^ 2)
  end;
  pack;
end;
x3arrAlias = linspace(x3min, x3max, numSteps * aliasFact);
y3arrAlias = linspace(y3min, y3max, numSteps * aliasFact);
[x3mat, y3mat] = meshgrid(x3arr, y3arr);
[x3matAlias, y3matAlias] = meshgrid(x3arrAlias, y3arrAlias);
finalLonerAlias = interp2(x3mat, y3mat, finalLoner, x3matAlias, y3matAlias, 'nearest');
contourf(x3matAlias, y3matAlias, transpose(finalLonerAlias), [1 2 3]);
colormap([1 0 0; 0 1 0; 0 0 1])
title(['Color (z) axis represents the number of the body that is on its own at the end of the simulation. ', num2str(numSteps), ' steps.']);

orient landscape
print -dpsc2 -Plp3