function  Orbit(N)
u1 = [1.2; 0; 0; -1.0493575];
Ti1 = 6.19216933;

u4 = [0.994; 0; 0; -2.031732629557];
Ti4 = 11.124;

u5 = [0.994; 0; 0; -2.113898796645];
Ti5 = 5.437;

options = odeset('RelTol',1e-10,'AbsTol',1e-10,'Stats','on');%,'Refine',1);
w1 = 1/82.45;
w2 = 0.012277471;
sprintf('ode23 for 1st closed orbit, Ti1= 6.19216933*4') 
[T1, Y1] = ode23(@equ927,linspace(0,4*Ti1,N), u1, options,w1);
sprintf('ode45 for 1st closed orbit, Ti1= 6.19216933*4')
[T2, Y2] = ode45(@equ927,linspace(0,4*Ti1,N), u1, options,w1);
sprintf('ode23 for 2nd closed orbit, Ti4 = 11.124*4')
[T3, Y3] = ode23(@equ927,linspace(0,4*Ti4,N), u4, options,w2);
sprintf('ode45 for 2nd closed orbit, Ti4 = 11.124*4')
[T4, Y4] = ode45(@equ927,linspace(0,4*Ti4,N), u4, options,w2);
sprintf('ode23 for 3rd closed orbit, Ti5 = 5.437*4')
[T5, Y5] = ode23(@equ927,linspace(0,4*Ti5,N), u5, options,w2);
sprintf('ode45 for 3rd closed orbit, Ti5 = 5.437*4')
[T6, Y6] = ode45(@equ927,linspace(0,4*Ti5,N), u5, options,w2);

function up = equ927(t,u,w)
%w = 1/82.45;
v = 1 - w;
%------------------------------------------
u1    = u(1);
u2    = u(2);
u3    = u(3);
u4    = u(4);

r1    = sqrt((u1+w)^2 + u3^2);
r2    = sqrt((u1-v)^2 + u3^2);

up    = zeros(4,1);

up(1) = u2;
up(2) = 2*u4 + u1 - v*(u1+w)/r1^3 - w*(u1-v)/r2^3;
up(3) = u4;
up(4) = -2*u2 + u3 - v*u3/r1^3 - w*u3/r2^3;

