Function Can't work

I'm having issues with making a function that I have received as an answer in the support forum of Mathworks (link: https://uk.mathworks.com/matlabcentral/answers/395183-creating-a-set-of-equally-spaced-cylinders-with-common-centre).

As soon as I copy and paste the function on the editor view, and I define points P1 and P2 as recommended, even after defining nVertex and R, the function appear to be working (no error message is shown), but at the end of the code, if I type mesh(RX,RY,RZ), then I the error message shown is that I have not defined RX,RY or RZ. Any ideas how I can make this work - where I should type my variables and which these variables are going to be so that the function successfully plots a mesh plot? The code is as follows and has been provided by Jan from Mathworks.

%these are the four variables that I define at the start of the code (of course without the '%' symbol):
%P1 = [0,0,0];

%P2 = [100,100,0];

%R = 2;

%nVertex = 20;

function H = DrawCylinder(P1, P2, R, nVertex, FaceColor, FaceLite)

[RX, RY, RZ] = CylinderCoor(P1, P2, R, nVertex);

H = surface(RX, RY, RZ, ...

    'FaceColor',        FaceColor, ...

    'FaceLighting',     FaceLite, ...

    'EdgeColor',        'none', ...

    'AmbientStrength',  0.5, ...

    'DiffuseStrength',  0.4, ...

    'BackFaceLighting', 'lit', ...

    'Clipping',         'off');

end

function [RX, RY, RZ] = CylinderCoor(P1, P2, R, nVertex)

smallVal = 1.49e-008;  % SQRT(EPS)

% Vector in the center of the cylinder:

nL12 = P2 - P1;

nL12 = nL12 ./ sqrt(sum(nL12 .* nL12, 2));

% Base vectors of a circle around the center with normal nL12:

cX    = [nL12(:, 2), -nL12(:, 1)];

LencX = sqrt(cX(:, 1) .* cX(:, 1) + cX(:, 2) .* cX(:, 2));

badInd        = (not(isfinite(LencX)) | LencX <= smallVal);

LencX(badInd) = 1.0;

cX(badInd, 1) = 1.0;

cX(badInd, 2) = 0.0;

cX       = cX ./ [LencX, LencX];

cX(:, 3) = 0;

% cY as normalized cross product of nL12 and cX:

cY = [-nL12(:, 3) .* cX(:, 2), ...

        nL12(:, 3)  .* cX(:, 1), ...

        nL12(:, 1)  .* cX(:, 2) - nL12(:, 2) .* cX(:, 1)];

cY = cY ./ sqrt(sum(cY .* cY, 2));

% Circle with nVertex points:

% theta = (0:nVertex) * (2 * pi / nVertex);

theta    = 0:(6.283185307179586 / nVertex):6.283185307179586;

costheta = R * cos(theta);

sintheta = R * sin(theta);

sintheta(nVertex + 1) = 0;  % Delete rounding errors (?!)

% Dyadic products: [M x 1] * [1 x N] = [M x N]

DX  = transpose(cX(:, 1) * costheta + cY(:, 1) * sintheta);

DY  = transpose(cX(:, 2) * costheta + cY(:, 2) * sintheta);

DZ  = transpose(cX(:, 3) * costheta + cY(:, 3) * sintheta);

P1T = transpose(P1);

P2T = transpose(P2);  

RX(2, :, :) = bsxfun(@plus, DX, P2T(1, :));

RY(2, :, :) = bsxfun(@plus, DY, P2T(2, :));

RZ(2, :, :) = bsxfun(@plus, DZ, P2T(3, :));

RX(1, :, :) = bsxfun(@plus, DX, P1T(1, :));

RY(1, :, :) = bsxfun(@plus, DY, P1T(2, :));

RZ(1, :, :) = bsxfun(@plus, DZ, P1T(3, :));

end

%mesh (RX,RY,RZ) (error message for undefined variables shown)