d1 = 1;
d2 = 10;
mu1 = 10;
mu2 = 1;
i = 0;
xvec = -20:0.1:20;
yvec = -20:0.1:20;
for x = xvec,
   i = i + 1;
   j = 0;
   for y = yvec,
      j = j  + 1;
      a = orbit(0,mu1,mu2,d1,d2,[x y 0],[0 0 0]);
      xa(i,j) = a(1);
      ya(i,j) = a(2);
   end
end
mat = 1000*sqrt(xa.^2+ya.^2);
idx = find(mat>64);
mat(idx) = 0;
image(yvec,xvec,mat');
s = jet(64);
s(1,:) = [1,1,1];
colormap(s); axis('equal');