function [a, t] = solsym2(x1,y1,z1,vx1,vy1,vz1,x2,y2,z2,vx2,vy2,vz2,t0,tf)
%input initial x and y position and their velocity components vx,vy
% and the start and finish times t0, tf

GM = 2.959122082856202e-04; % GMsun in AU^3/day^2
% GM = (1.32712440018e20 (m3s-2)/(1.49597870691e11(m/AU))^3)*86400(sec/day)^2;

% planet masses
m = [GM 0.000002*GM];
X = [x1,y1,z1,vx1,vy1,vz1,x2,y2,z2,vx2,vy2,vz2];
tspan = [t0,tf];

options = odeset('Stats','on', 'RelTol', 1e-15 );

% use ode45 function to solve differential equations
[t,a] = ode45(@grav2,tspan,X, options, m);

whos a

% plot data in x,y plane
plot3(a(:,1),a(:,2),a(:,3),a(:,7),a(:,8),a(:,9));
