clc;
 clear all;
 z = -25:0.1:400;
 n = size(z);
 for i = 1:n(2)
     if     (z(i) > 0)
         S(i) = (z(i)^0.5 -sin(z(i)^0.5))/z(i)^1.5;
         C(i) = (1 - cos(z(i)^0.5))/z(i);
     elseif (z(i) < 0)
         S(i) = (sinh((-z(i))^0.5) - (-z(i))^0.5)/((-z(i))^1.5);
         C(i) = (cosh((-z(i))^0.5) - 1)/(-z(i));
     else
         S(i) = 1/6;
         C(i) = 0.5;
     end
 end
% Plot
figure(1);
scatter(z(1:500),C(1:500),'.b');
hold on;grid on;
scatter(z(1:500),S(1:500),'.r');
xlabel('z');
ylabel('C(z),S(z)');
legend('C(z)','S(z)');
title('Stumpff functions C(z),S(z)');
text(10,0.4,'smallsats.org','Color',[0 0 0], 'VerticalAlignment','middle',...
'HorizontalAlignment','left','FontSize',14 );
hold off;

figure(2);
scatter(z(501:n(2)),C(501:n(2)),'.b');
hold on;grid on;
scatter(z(501:n(2)),S(501:n(2)),'.r');
xlabel('z');
ylabel('C(z),S(z)');
legend('C(z)','S(z)');
title('Stumpff functions C(z),S(z)');
text(250,0.0010,'smallsats.org','Color',[0 0 0], 'VerticalAlignment','middle',...
'HorizontalAlignment','left','FontSize',14 );
hold off;

Equations from the book Orbital Mechanics for Engineering Students, Second Edition, Aerospace Engineering