How can I add a trailing singleton dimension to a matrix

You can use permute to inject a trailing default singleton dimension into a non-trailing position. For instance:

A = rand(50,50); % produces a 50 x 50 matrix
B = permute(A, [1 3 2]); % produces a 50 x 1 x 50 array

Here, 3 refers to the third dimension in the array, an implicit dimension of size 1.

The trailing singleton is dropped by MATLAB even from the point of view of the MEX API. The docs for mxSetDimensions say as much:

MATLAB® automatically removes any trailing singleton dimensions specified in the dims argument. For example, if ndim equals 5 and dims equals [4 1 7 1 1], the resulting array has the dimensions 4-by-1-by-7.

Again, that's from the mxSetDimensions docs.

Here's a test to try this with mxSetDimensions or mxCreateNumericArray:

// mexSizeTest.cpp
#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[])
{
    if (nrhs!=1)
        mexErrMsgTxt("Just one input please.");

    mxArray *M = mxDuplicateArray(prhs[0]);
    const mwSize ndimsOrig = mxGetNumberOfDimensions(M);
    const mwSize *dimsOrig = mxGetDimensions(M);

    mexPrintf("Number of dimensions (input): %d\n\n", ndimsOrig);

    const mwSize ndims = ndimsOrig + 1;
    mwSize *dims = (mwSize*) mxMalloc(ndims*sizeof(mwSize));
    for (int i=0; i<ndimsOrig; ++i) dims[i] = dimsOrig[i];
    dims[ndims-1] = 1;

    mexPrintf("Target dimensions: [");
    for (int i=0; i<ndims-1; ++i) mexPrintf("%d ",dims[i]);
    mexPrintf("%d]\n\n",dims[ndims-1]);

    mexPrintf("Reshaping to #dims = %d with trailing singleton.\n", ndims);
    mxSetDimensions(M, dims, ndims);
    mexPrintf("Number of Dimensions: %d\n\n", mxGetNumberOfDimensions(M));

    // Let's be dangerous to see if the 1 is just hiding under the hood
    const mwSize *dimsSet = mxGetDimensions(M);
    mexPrintf("Being dangerous: %d\n\n", dimsSet[ndims-1]); // !!!
    mxDestroyArray(M);

    mexPrintf("Creating fresh mxArray of #dims = %d with trailing singleton.\n",
            ndims);
    M = mxCreateNumericArray(ndims, dims, mxDOUBLE_CLASS, mxREAL);
    mexPrintf("Number of Dimensions: %d\n",mxGetNumberOfDimensions(M));
    mxDestroyArray(M); mxFree(dims);
}

The MATLAB test:

>> M = rand(24,55,1);
>> size(M)
ans =
    24    55
>> ndims(M)
ans =
     2
>> size(M,454235) % note we can check the 454235th dimension
ans =
     1

Side note with that size(M,454235) test: This is what the ndim docs mean when they say Trailing singleton dimensions are ignored. They're really not ignored, they just aren't really there!

The MEX test (mexSizeTest.cpp):

>> mexSizeTest(M)
Number of dimensions (input): 2

Target dimensions: [24 55 1]

Reshaping to #dims = 3 with trailing singleton.
Number of Dimensions: 2

Being dangerous: -994713024

Creating fresh mxArray of #dims = 3 with trailing singleton.
Number of Dimensions: 2

I suggest to adapt your code to handle the case where mxGetNumberOfDimensions returns 2.