Many thanks for the changes!
I'm ok with adapting the description and also with keeping lines 69 - 71 active.

In the meantime I detected one more smaller issue:
At the moment, pammod will change the orientation between input and output vector, esp. when input is a column vector. Adding the following after line 71 would also solve this:

y = reshape(y,size(x));

You are right, either using the complex function ore taking a very small phi would produce the desired complex output.
So, I'll keep it up to you to decide if pammod would be changed.

When the current implementation will be kept, lines 69 to 71 would never be invoked and, hence, could be disabled.

I fully agree, that from conceptual point of view the current implementation delivers correct results when looking purely at real value modulation. However, the problem occurs when putting it together with the awgn function. In this case pure real values will generate the noise only on the real axis. This results in the lines on the x-axis visible in scatterplot. When having complex values (even with imaginary part equal to zero as provided by Matlab), awgn will generate the noise in both domains. This is visible as the circles around the constellation points in the scatterplot. With the same dB, the real only values produce more errors than the complex values. This behaviour is in principle also ok, but would now result in a difference between Octave and Matlab. To illustrate this, I attach another file, which produces a BER-plot (please note that qamdemod is not working correctly in current Octave implementation). With this, Matlab is producing the values, which are quite aligned with the theoretical model, while Octave differs by 3 dB. Interestingly, when changing phi by a very small amount (e.g. setting to eps), the results will fully align to the Matlab case.

(file #56720)

When simulating a simple communication the current behaviour produces together with the awgn function too many errors (cf. the enclosed test file). In fact, the SNR would need to be shifted by 3 dB. Only with complex values the correct number of errors are produced that fit to the theoretical behaviour of PAM.
BTW, my Matlab (at home I'm using the 2017 version, but I believe that this has not changed), produces with the same file complex modulation symbols.

(file #56717)

As in the associated Matlab function, pammod should provide complex output even in case of pure real values. An implementation example is given in the attached file (starting at line 69).