function null_img
% Schwarzchild time-like geodesic
% with imaginary eccentricities
e = .01;     % eccentricity 
M = 1/3;    % mass of the black hole
l = 1;   % latus rectum

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

options = odeset('Refine',38);

xispan = [2*atan(-1/e) pi];    %bound xi_inf = 
                               %2*atan(-(6*mu-1)/2/mu/e)
[fi, xi] = ode45(@fi, xispan, 0, options, e, ...
    mu);

polar(fi, l./(1+e*tan(xi/2)))
%----------------------------------------------

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