function slpotent
% Reissner-Nordstrom space-like geodesic
% plotting of potential
% M is a mass of black hole
Q = 0.8;                 % charge of black hole

r = .005:0.005:10;

subplot 311
M = 10/14;
plot(r, potentia(M, Q, 3, r),  r, potentia(M, Q, 5, r), 'b', ...
    r, potentia(M, Q, 10, r), 'b', r, potentia(M, Q, 15, r), ...
    'b', 'LineWidth', 2)
axis([0 10 -10 40])
set(gca, 'Box', 'off')
title('Case when {\it M} = 10/14 <{\it M}_{cr}', 'Fontsize', 15)
xlabel('{\it r}', 'FontSize', 16, 'Position', [9.5 -11])
ylabel('{\it V}', 'FontSize', 16, ...
    'Position', [-1 27], 'Rotation', 0)
GTEXT('a', 'FontWeight', 'bold')
GTEXT('b', 'FontWeight', 'bold')
GTEXT('c', 'FontWeight', 'bold')
GTEXT('d', 'FontWeight', 'bold')
subplot 312;
plot(r, potentia(mass(3, Q), Q, 3, r), 'b', r, ...
    potentia(mass(5, Q), Q, 5, r), 'b', r, potentia(mass(10, ...
    Q), Q, 10, r), 'b', r, potentia(mass(15, Q), Q, 15, r), ...
    'b', 'LineWidth', 2)
axis([0 10 -5 40])
set(gca, 'Box', 'off')
title('Special case when {\it M} ={\it M}_{cr}', ...
    'Fontsize', 15)
xlabel('{\it r}', 'FontSize', 16, 'Position', [9.5 -6])
ylabel('{\it V}', 'FontSize', 16, ...
    'Position', [-1 30], 'Rotation', 0)
GTEXT('a', 'FontWeight', 'bold')
GTEXT('b', 'FontWeight', 'bold')
GTEXT('c', 'FontWeight', 'bold')
GTEXT('d', 'FontWeight', 'bold')
subplot 313
M = 12/14;
plot(r, potentia(M, Q, 3, r), 'b', r, potentia(M, Q, 5, r), ...
    'b', r, potentia(M, Q, 10, r), 'b', r, ...
    potentia(M, Q, 15, r), 'b', 'LineWidth', 2)
axis([0 10 -70 40])
set(gca, 'Box', 'off')
title('Case when {\it M} = 12/14 >{\it M}_{cr}', 'Fontsize', 15)
xlabel('{\it r}', 'FontSize', 16, 'Position', [9.5 -70])
ylabel('{\it V}', 'FontSize', 16, ...
    'Position', [-1 3.8], 'Rotation', 0)
GTEXT('a', 'FontWeight', 'bold')
GTEXT('b', 'FontWeight', 'bold')
GTEXT('c', 'FontWeight', 'bold')
GTEXT('d', 'FontWeight', 'bold')

%-------------------------------------------------------------
function [V] = potentia(M, Q, L, r)
% Potential V for parametres M, L and Q (L is angular momentum
% about an axis normal to the invariant plane, M is mass,
% Q is charge)
V = (r.^2-2*M*r+Q^2) ./r.^2 .*(L^2./r.^2-1);

function [M] = mass(L, Q)
% Mass for special case, when f(u) has an inflexion
M = sqrt(6*(-Q^4+14*Q^2*L^2-L^4) ...
    +2*sqrt(3*(3*L^2-Q^2)*(L^2+5*Q^2)^3)) /12/ L;