function spacelik
% Schwarzchild  geodesic
% of the first and the second kind
e = 3/2;       % eccentricity 5 - space-like  
               % 1/2 - bound, 3/2 unbound - time-like
M = 2/14;      % mass of the black hole
l = 4.5;       % latus rectum 4.5 - unbound,   3 - bound

mu = M/l;
%1 -mu*(3+e^2)    % > 0
%l/(e+1)
k = 2*sqrt(mu*e/(1-6*mu+2*mu*e));

chispan = [0 acos(-1/e)]; % [0 acos(-1/e)] - unbound   
                          % [0 12*pi]  - bound 
options = odeset('Refine',38);
[fi1, chi] = ode45(@fi1, chispan, -1, options, e, ...
    mu, k); % chi0 = 0 - unbound,  0 - bound

xispan = [0 pi];
[fi2, xi] = ode45(@fi2, xispan, 0, options, e, ...
    mu, k);

polar(fi1, radius1(chi, l, e))
hold on
%polar(-fi1, radius1(chi, l, e))     % <--- for unbound
%hold on
polar(fi2, radius2(xi, l, e, M))
title(['M = ' num2str(M) ', l = ' num2str(l), ...
        ', e = ' num2str(e)], 'Fontsize' , 15);
%----------------------------------------------

function f1 = fi1(fi, chi, e, mu, k)
% ODE for phi coordinate for the geodesics of the first kind
% f1 is d chi / d phi
f1 = sqrt(1-6*mu+2*mu*e) *sqrt(1-k^2*cos(chi/2)^2);

function f2 = fi2(fi, xi, e, mu, k)
% ODE for phi coordinate for the geodesics of the second kind
% f2 is d xi / d phi
f2 = sqrt(1-6*mu+2*mu*e) *sqrt(1-k^2*sin(xi/2)^2);