eta = 2*asinh(-1):0.05:2*asinh(1);
r_i = 1; %r_i = M = 1, (E = 1) initial conditions
u = sqrt(1-radius(eta, r_i)/2) .*(exp(radius(eta, r_i)/4)) ...
    .*sinh(time(eta, r_i)/4);
v = sqrt(1-radius(eta, r_i)/2) .*(exp(radius(eta, r_i)/4)) ...
    .*cosh(time(eta, r_i)/4);
a = [-2 2]; %for plotting r = 2M

% singularity r = 0:
u2 = sinh(-pi/2:0.05:pi/2);
v2 = cosh(-pi/2:0.05:pi/2);

plot(a, a, ':', a, -a, ':', u, v)
hold on
plot(u2, v2, -u2, -v2, 'LineWidth', 2)
axis([-2 2 -2 2])
set(gca, 'Box', 'off')
title('Kruskal-Szekeres coordinates', 'Fontsize', 15)
xlabel('u', 'FontAngle', 'italic', 'FontSize', 16, ...
    'Position', [1.98 -2.2])
ylabel('v', 'FontAngle', 'italic', 'FontSize', 16, ...
    'Position', [-2.2 1.9], 'Rotation', 0)
x = [-2 2];
y = [0 0];
line(x, y)
line(y, x)
text(.7,.9, '{\it r} = 2{\it M}, {\it t} = +\infty', ...
    'Rotation', 38, 'VerticalAlignment', 'Cap', ...
    'FontWeight', 'bold')
text(-1.3,1.5, '{\it r} = 2{\it M}, {\it t} = -\infty', ...
    'Rotation', -38, 'VerticalAlignment', 'Cap', ...
    'FontWeight', 'bold')
text(-.3, 1.5, '{\it r} = 0', 'FontWeight', 'bold')
text(-.3, -1.5, '{\it r} = 0', 'FontWeight', 'bold')