function spacelik
% Kerr spacelike non-equatorial
% unstable circular geodesics
M = 1;      % mass of the black hole
a = .96;    % angular momentum per unit mass of the black hole
r_s = 1.85; % radius of the circular orbit 
E = 1;      % energy of the particle

Delta_s = r_s^2 -2*M*r_s +a^2;  % horizon function for circle
alpha = a *sqrt(1+1/E^2);
eta_s = (2*r_s^3*M*Delta_s/(r_s-M) ...
    *(1-sqrt(1+1/E^2*(1-M/r_s))) +r_s^4*(r_s-M) ...
    -r_s^2/E^2*(r_s*(r_s-2*M)^2-a^2*M)) /a^2 /(r_s-M);
xi_s = (M*(r_s^2-a^2)+r_s*Delta_s*sqrt(1+1/E^2*(1-M/r_s))) ...
    /a /(r_s-M);
acos((sqrt((sqrt((xi_s^2+eta_s-alpha^2)^2+4*alpha^2*eta_s))/2/alpha^2)))
eta_s-(xi_s
tauspan = [.1  3];
initcond = [2 1.5 0];
options = odeset('Refine', 38, 'Mass', @Mass, ...
    'MStateDependence', 'strong', 'MvPattern', MvPat, ...
    'MassSingular', 'no');

[tau, y] = ode15s(@nullKerr, tauspan, initcond, options, ...
    a, xi_s, eta_s, M, E, r_s, alpha);

subplot(2,1,1)
plot(y(:,1).*cos(y(:,3)).*sin(y(:,2)), ...
   y(:,1).*cos(y(:,2)))
title('Projection on the (r, \theta) plane', ...
    'Fontsize', 15);

subplot(2,1,2)
plot(y(:,1).*cos(y(:,3)).*sin(y(:,2)), ...
    y(:,1).*sin(y(:,3)).*sin(y(:,2)))

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(['Space-like circular orbit with (r_s = ', ...
    num2str(r_s), ')'], 'Fontsize', 15);
%-------------------------------------------------------------

function dy = nullKerr(tau, y, a, xi_s, eta_s, M, E, r_s, ...
    alpha)
% 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) = (y(1)-r_s)*sqrt((1+1/E^2)*y(1)^2 ...
    +2*r_s*y(1)*(1+1/E^2*(1-M/r_s)) -a^2*eta_s/r_s^2);
dy(2) = sqrt(eta_s +alpha^2*cos(y(2))^2 ...
    -xi_s^2*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_s*csc(y(2))^2;

function Mass = Mass(tau, y, a, xi_s, eta_s, M, E, r_s, alpha)
% 'mass' matrix for ODE solver r^2 + a^2 cos theta
rho = (1/y(1)+r_s)^2 +a^2*y(2);
Delta = (1/y(1)+r_s)^2 -2*M*(1/y(1)+r_s) +a^2;
rho_vek = [rho^2 rho^2 rho^2*Delta];
Mass = diag(rho_vek);

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