function imaginar
% Schwarzchild geodesic
% with the imaginary eccentrities
e = 1;     % eccentricity
M = 2/14;  % mass of the black hole: 2/14 -- space-like
                                    
l = 1;   % latus rectum: 0.16 -- space-like
                            % 1 -- time-like
mu = M/l;                            % < .25
% Delta = sqrt((6*mu-1)^2+4*mu^2*e^2);
% k = sqrt((Delta +6*mu-1)/2/Delta) % > 1
% 1-3*mu+mu*e^2                     % > 0
% 2*ellipke(k)/sqrt(Delta)
% 2*atan(-(6*mu-1)/2/mu/e)
% 1/(1/2/M-2/l)                         % > 0
 
xispan = [1 5.9];
                    % [-1.57 .1] -- space-like
                    % [1 5.9] -- time-like
options = odeset('Refine',38);

[fi, xi] = ode45(@fi, xispan, 1, ...
    options, e, ...    
    mu);        % init. cond: -1 -- space-like, 1 -- time-like

polar(fi, i_radius(xi, l, e))
title(['M = ' num2str(M) ', l = ' num2str(l) ...
        ', e = ' num2str(e) '\cdot i'], 'Fontsize', 15);
%----------------------------------------------

function f = fi(fi, xi, e, mu)
% ODE for phi coordinate for the geodesics
% with imaginary eccentricity
% f1 is d xi / d phi
f = sqrt(2*(6*mu -1 +2*mu*e*sin(xi) +(6*mu-1)*cos(xi)));