function null_eq
% Kerr null geodesic
% in equatorial plane
% flyby orbits
a = .96;     % angular momentum per unit mass 
M = 1;   % mass of the black hole
D = 2;   % impact parameter

options = odeset('Refine',38);

uspan1 = [.33 0];
[fi1, u1] = ode45(@fi, uspan1, pi/3, options, a, ...
    M, D);
uspan2 = [1 0];
[fi2, u2] = ode45(@fi, uspan2, pi/3, options, a, ...
    M, D);
%plot(fi1, u1, fi2, u2)
polar(fi2, 1./u2)
hold on
polar(fi1, 1./u1)
hold on
polar(-fi1, 1./u1)

%----------------------------------------------

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