How to create a surface plot of a function that uses colon operators?

I have a following function:

function p =  dist(n, theta, M, Ns, Nb, G, kappa)
    Ntheta =  55;
    logterm1 = sum( log(1: n + M - 1) ) - sum( log (1: M - 1)) - sum( log (1:n));
    logterm2 =  n.*(log (Ntheta) - log(1 + Ntheta));
    logterm3 = M*log(1 + Ntheta);
    total = logterm1 + logterm2 - logterm3;
    p = exp(total);
end


function PE = Error(x, y)

    p0 = x;
    p1 = 1-x;
    Ns = 1e-2;
    Nb =10;
    kappa = 1e-2;
    G = 1.01;
    M = 1e3;
    Ntheta0 =  150;
    Ntheta_pi =  33;
    sigma0 = sqrt( Ntheta0.*( Ntheta0 + 1));
    sigma_pi = sqrt( Ntheta_pi.*( Ntheta_pi + 1));
    
    narray = 1:floor(y);

    PX0LessthanTh = 0;
    for na = narray
        pp = dist(na, 0, M, Ns, Nb, G, kappa);
        PX0LessthanTh = PX0LessthanTh + pp;
    end

    PXpiLessthanTh = 0;
    for na = narray
        pp = dist(na, pi, M, Ns, Nb, G, kappa);
        PXpiLessthanTh = PXpiLessthanTh + pp;
    end
    PE = p0.*PX0LessthanTh + p1.*(1 - PXpiLessthanTh);

end

I am looking to create a surface plot of Error with respect to x and y. However, I because I am using a statement narray = 1:floor(y); I am getting an error:

 The following error was reported evaluating the function in FunctionLine
 update: Colon operands must be real scalars.

How can I modify my Error function above so that I can use 3D functions like surf.