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

%mu = M/l;

options = odeset('Refine',38);

uspan = [0 100];
[fi, u] = ode45(@f, uspan, 100, options, e, ...
    M, l);

polar(fi, l./u)
%----------------------------------------------

function f = f(fi, u, e, M, l)
% ODE for phi coordinate
% f is d u / d phi
f = sqrt(2*M) .*sqrt(u-1/2/M-2/l) ...
    .*sqrt((u-1/l).^2+e^2/l^2);