function plottime
% Schwarzchild space-like geodesic
% of the first and the second kind
e = 6;     % eccentricity 
M = 2/14;  % mass of the black hole
l = 4.5;   % latus rectum

mu = M/l;
k = 2*sqrt(mu*e/(1-6*mu+2*mu*e));

zac = -4.2 ;
timespan = [zac 0];  
chi0 = 0;
options = [];
[pt1, pchi] = ode45(@tau1, timespan, chi0, options, e, mu, ...
    k, l);
[st1, schi] = ode45(@t1, timespan, chi0, options, e, mu, k, l);

tauspan = [zac -3.937];   
xi0 = 0;
[pt2, pxi] = ode45(@tau2, tauspan, xi0, options, e, mu, k, l);
[st2, sxi] = ode45(@t2, timespan, xi0, options, e, mu, k, l);

subplot(211)
plot(pt1, radius1(pchi, l, e), pt2, radius2(pxi, l, e, M), ...
    'b', 'LineWidth', 2);
xlabel('\tau', 'FontAngle', 'italic','FontSize', 16, ...
    'Position', [-2.02 -0.4]);
ylabel('r', 'FontAngle', 'italic','FontSize', 16, ...
    'Position', [-4.61 2.8], 'Rotation', 0);
title('Proper time', 'Fontsize', 15);
text(-3.55, .6, 'first kind', 'FontWeight', 'bold');
text(-3.55, .35, 'second kind', 'FontWeight', 'bold');
axis([-4.5 -2 0 3]);
set(gca, 'Box', 'off');
subplot(212)
plot(st1, radius1(schi, l, e), st2, radius2(sxi, l, e, M), ...
    'b', 'LineWidth', 2 );
xlabel('t', 'FontAngle', 'italic','FontSize', 16, ...
    'Position', [-0.52 -0.4]);
ylabel('r', 'FontAngle', 'italic','FontSize', 16, ...
    'Position', [-4.7 2.8], 'Rotation', 0);
title('Schwarzschild time', 'Fontsize', 15);
text(-1.6, 1, 'first kind', 'FontWeight', 'bold');
text(-1.6, .75, 'second kind', 'FontWeight', 'bold');
set(line([-4.5 0], [2*M 2*M]), 'LineStyle', ':');
axis([-4.5 -0.5 0 3]);
set(gca, 'Box', 'off');
%----------------------------------------------

function f1 = tau1(tau, chi, e, mu, k, l)
% ODE for proper time for the geodesics of the first kind
% f1 is d chi / d tau
f1 = 2*mu/k/l *sqrt(e/(3*mu+mu*e^2-1)) *(1+e*cos(chi))^2 ...
    *sqrt(1-k^2*cos(chi/2)^2);

function g1 = t1(t, chi, e, mu, k, l)
% ODE for schwarzschild time for the geodesic of the first kind
% g1 is d chi / d t
g1 = 2*mu*sqrt(e) /k/l *(1-2*mu-2*mu*e*cos(chi)) ...
    *(1+e*cos(chi))^2 /sqrt(1-4*mu+4*mu^2-4*mu^2*e^2) ...
    *sqrt(1-k^2*cos(chi/2)^2);

function f2 = tau2(tau, xi, e, mu, k, l)
% ODE for proper time for the geodesics of the second kind
% f2 is d xi / d tau
f2 = sqrt(e) /2/mu/k/l ...
    *(1+4*mu+(1-6*mu-2*mu*e)*tan(xi/2)^2)^2 ...
    *sqrt(3*mu+mu*e^2-1) *sqrt(1-k^2*sin(xi/2)^2);

function g2 = t2(t, xi, e, mu, k, l)
% ODE for schwarzschild time for the geodesic of the second kind
g2 = sqrt(e-e*k^2*sin(xi/2)^2) /2/mu/k/l ...
    /sqrt(1-4*mu+4*mu^2-4*mu^2*e^2) ...
    *(1-4*mu+(1-6*mu-2*mu*e)*tan(xi/2)^2)^2 ...
    *(4*mu -(1-6*mu-2*mu*e)*tan(xi/2)^2);