% initial positions
ip1 = [1 0]; ip2 = [-0.5 0.86]; ip3 = [-0.5 -0.86];
% initial velocities
iv1 = [-1 0]; iv2 = [-0.7 -0.7]; iv3 = [-0.7 0.7];
% input for initial righthandside vector
initz = [ip1(1) iv1(1) ip1(2) iv1(2) ip2(1) iv2(1) ip2(2) iv2(2) ip3(1) iv3(1) ip3(2) iv3(2)]
% computing the initial righthandside vector 
r = planar3rhsOctave(initz,0)
% solving the IVP :
tspan = [0:0.03:1];
z = lsode('planar3rhsOctave',initz,tspan);
% extracting the trajectories
x1 = z(:,1); y1 = z(:,3);
x2 = z(:,5); y2 = z(:,7);
x3 = z(:,9); y3 = z(:,11);
% plotting the trajectories
hold on
plot(x1,y1,'r')
plot(x2,y2,'g')
plot(x3,y3,'b')