function hkr
% zkouska clanku hkr

m = 1;  % mass of the black hole
l = .75;   % angular momentum
gama = sqrt(.36); %energie

options = [];
fispan = [0 pi/3.19];
[fi1, r] = ode45(@eq10, fispan, 5, options, l, ...
     gama, m);
[fi2, r2] = ode45(@eq10b, [0 1.4*pi], 0.1, options, l, ...
     gama, m);

% r=0.5;
%gama^2-2*m/r-l^2/r^2+1+2*m*l^2/r^3
%polar(fi1, r)
%hold on
%polar(-fi2, r2)
%hold on
%polar(-fi1, r)
%hold on
polar(-fi2, r2)
%----------------------------------------------

function drdfi = eq10(fi1, r, l, gama, m)
% ODE for phi and r coordinates for the geodesics of the first kind
% drdfi is d r / d phi
drdfi = -r.^2/l .*sqrt(gama^2-2*m./r-l^2/r.^2+1+2*m*l^2./r.^3);

function drdfi = eq10b(fi2, r, l, gama, m)
% ODE for phi and r coordinates for the geodesics of the first kind
% drdfi is d r / d phi
drdfi = r.^2/l .*sqrt(gama^2-2*m./r-l^2./r^2+1+2*m*l^2./r^3);