We start with some of the simplest possible NAG calls. Among many other special functions, the S chapter contains
the Bessel functions Y0(x), Y1(x), J0(x) and J1(x). Following the usual NAG naming convention, these are in routines
s17ac, s17ad, s17ae and s17af respectively.
Here's the code, including some plotting commands:

% Evaluate the Bessel functions at a set of points in x
x = [0.125, 0.25 : 0.25 : 40];
for i = 1 : length(x)
[y0(i), ifail] = s17ac(x(i));
[y1(i), ifail] = s17ad(x(i));
[j0(i), ifail] = s17ae(x(i));
[j1(i), ifail] = s17af(x(i));
end
% Plot the results. Note that we omit the first two values of Y1(x)
% because they are large and would spoil the scaling of the graph.
hold on;
plot(x, y0, 'Color', 'red');
plot(x(3:length(x)), y1(3:length(x)), 'Color', 'green');
plot(x, j0, 'Color', 'blue');
plot(x, j1, 'Color', [1.0 1.0 0.0]);


and the results are shown in this picture:
Figure: the Bessel functions