function pl_itime
% Schwarzchild space-like geodesic
% with immaginary eccentricity
e = 1;     % eccentricity 
M = 2/14;  % mass of the black hole
l = 0.16;   % latus rectum

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

zac = -4.1;
timespan = [zac 1.055];
xi0 = -1;
options = [];
[pt, pxi] = ode45(@tau, timespan, xi0, options, e, mu, k, l);
[st, sxi] = ode45(@t, timespan, xi0, options, e, mu, k, l);

subplot(211)
plot(pt, i_radius(pxi, l, e), 'LineWidth', 2);
xlabel('\tau', 'FontAngle', 'italic','FontSize', 16, ...
    'Position', [1.5 -0.05]);
ylabel('r', 'FontAngle', 'italic','FontSize', 16, ...
    'Position', [-5.1 .35], 'Rotation', 0);
title('Proper time', 'Fontsize', 15);
axis([-4.5 2 0 0.4]);
set(gca, 'Box', 'off');
subplot(212)
plot(st, i_radius(sxi, l, e), 'LineWidth', 2);
xlabel('t', 'FontAngle', 'italic','FontSize', 16, ...
    'Position', [1.5 0.27]);
ylabel('r', 'FontAngle', 'italic','FontSize', 16, ...
    'Position', [-5.1 .37], 'Rotation', 0);
title('Schwarzschild time', 'Fontsize', 15);
set(line([-4.5 0], [2*M 2*M]), 'LineStyle', ':');
axis([-4.5 2 0.28 .38]);
set(gca, 'Box', 'off');
%----------------------------------------------

function f = tau(tau, xi, e, mu, k, l)
% ODE for proper time for the geodesics with immaginary
% eccentricity
% f is d xi / d tau
f = mu *sqrt(2) /(mu*(3-e^2)-1) *(1+e*tan(xi/2))^2 ...
    /sqrt(6*mu -1 +2*mu*e*sin(xi) +(6*mu-1)*cos(xi));

function g = t(t, xi, e, mu, k, l)
% ODE for schwarzschild time for the geodesic with immaginary
% eccentricity
% g is d xi / d t
g = (1+e*tan(xi/2))^2/l *(1 -2*mu*(1+e*tan(xi/2))) ...
    *sqrt(2*mu/(4*mu^2*e^2+4*mu^2-4*mu+1)) ...
    *sqrt(6*mu -1 +2*mu*e*sin(xi) +(6*mu-1)*cos(xi));