function tlpotent
% Reissner-Nordstrom time-like geodesic
% plotting of potential
% M is a mass of black hole
Q = .8;                 % charge of black hole

r = .005:0.005:600;

%subplot 311
M = 10/14;
%plot(r, potentia(M, Q, 3, r),  r, potentia(M, Q, 5, r), r, ...
%    potentia(M, Q, 10, r), ...
%    r, potentia(M, Q, 15, r), 'LineWidth', 2)
%axis([0 600 0.94 1.02])
%set(gca, 'Box', 'off')
%title('Case when {\it M} = 10/14 <{\it M}_{cr}', ...
%    'Fontsize', 15)
%xlabel('{\it r}', 'FontSize', 16, 'Position', [580 .925])
%ylabel('{\it V}', 'FontSize', 16, ...
%    'Position', [-60 1.01], 'Rotation', 0)
%GTEXT('a', 'FontWeight', 'bold')
%GTEXT('b', 'FontWeight', 'bold')
%GTEXT('c', 'FontWeight', 'bold')
%GTEXT('d', 'FontWeight', 'bold')
%subplot 312;
plot(r, potentia(mass_p(.3, Q), Q, .3, r), r, ...
    potentia(mass_p(.5, Q), Q, .5, r), r, potentia(mass_p(1, ...
    Q), Q, 1, r), r, potentia(mass_p(1.5, Q), Q, 1.5, r),...
    'LineWidth', 2)
axis([0 10 -5 4])
%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), r, potentia(M, Q, .5, r), r, ...
%    potentia(M, Q, 1.0, r), r, ...
%    potentia(M, Q, 5, r), 'LineWidth', 2)
%axis([0. 2 -5 5])
%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_p(L, Q)
% Mass (biger) 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;