function nullneta
% Kerr null non-equatorial geodesics
% with eta < 0
eta = -0.3;     % related to Carter constant, negative
xi = -0.05;     % related to particle's angular momentum
M = 1;          % mass of the black hole
a = .8;         % angular momentum per unit mass of black hole

tauspan = [-0.3  -.85101925];
initcond = [1 0.6265 1];
options = odeset('Refine', 38, 'Mass', @Mass, ...
    'MStateDependence', 'strong', 'MvPattern', MvPat, ...
    'MassSingular', 'no');

[tau, y] = ode15s(@nullKerr, tauspan, initcond, options, ...
    a, xi, eta, M);

plot(y(:,1).*cos(y(:,3)), y(:,1).*sin(y(:,3)))
title('Projection on the (r, \phi) plane', ...
    'Fontsize', 15);
figure
plot3(y(:,1).*cos(y(:,3)).*sin(y(:,2)), ...
    y(:,1).*sin(y(:,3)).*sin(y(:,2)), y(:,1).*cos(y(:,2)))
set(findobj('Type', 'line'), 'LineWidth', 2)
title(['Null orbit with negative Carter constant (\eta = ', ...
    num2str(eta), ')'], 'Fontsize', 15);
%-------------------------------------------------------------

function dy = nullKerr(tau, y, a, xi, eta, M)
% dy(1) is d r / d tau
% dy(2) is d theta / d tau
% dy(3) is d phi / d tau
dy = zeros(3, 1);
dy(1) = sqrt(y(1)^4 +(a^2-xi^2-eta)*y(1)^2 ...
    +2*M*(eta+(xi-a)^2)*y(1) -a^2*eta);
dy(2) = sqrt(eta +(a*cos(y(2)))^2 -(xi*cot(y(2)))^2);
dy(3) = 2*a*M*y(1) ...
    +((y(1)^2+a^2*cos(y(2)))^2-2*M*y(1))*eta*csc(y(2))^2;

function Mass = Mass(tau, y, a, xi, eta, M)
% 'mass' matrix for ODE solver r^2 + a^2 cos theta
rho = y(1)^2 +a^2*cos(y(2));
Delta = y(1)^2 -2*M*y(1) +a^2;
rho_vec = [rho^2 rho^2 rho^2*Delta];
Mass = diag(rho_vec);

function S = MvPat
% sparsity pattern
% S(i, j) = 1 if for any k, the (i, k) component of 
% Mass(tau, y) depends on component j of y, and 0 otherwise
S = sparse(3, 3);
for i = 1:2
    for j = 1:3
        S(i, j) = 1;
    end
end
S(1, 1) = 1;