Given odd number $n$ count the bases to which $n$ is Euler pseudoprime

Question

As the title says we are given an odd number $n$ and wish to find the number of bases $b$ such that $n$ is an Euler pseudoprime; That is, $\gcd(b,n)=1$ and $b^{(n-1)/2} \equiv \left( \frac{b}{n} \right) \mod n$. I wrote a simple function in Matlab to compute the Jacobi symbol and used this in a simple for-loop over the numbers that are coprime to $n$. I tested it for some random numbers and although the numbers seem to check out I feel that it somehow gives me the wrong answer. Here is my function:

function [count] = EulerPP(n)
count=0;
for i=2:(n-1)
    if gcd(i,n)==1
        if mod(modexp(i,(n-1)/2,n)-JacSym(i,n),n)==0
            count=count+1;
        end
    end
end
end

Here is the input/output for some numbers:

$\begin{align*} 1387 &\mapsto 126 \\ 1189 &\mapsto 6 \\ 1537 &\mapsto 7 \\ 105 &\mapsto 7 \\ 781 &\mapsto 55 \\ \end{align*} $

Also, is it possible to calculate this in a simpler way (i.e. by not checking all numbers)?

Edit: I have also included the code for the functions modexp and JacSym:

function [s] = JacSym(m,n)
%Jacobi symbol of (m/n) where n is always odd.

m=mod(m,n);

if mod(n,2)==0
    s=0;
elseif m==1;
    s=1;
elseif m==0
    s=0;

elseif mod(m,2)==0
    if abs(mod(n,8))==1
        s=JacSym(m/2,n);
    else
        s=-JacSym(m/2,n);
    end
else
    if mod(n,4)==3 && mod(m,4)==3
        s=-JacSym(n,m);
    else
        s=JacSym(n,m);
    end
end    
end


function result = modexp (x, y, n)
    %anything raised to 0th power = 1 so return 1
    if (y == 0)
        result = 1;
        return;
    end

    %recurse
    z = modexp(x, floor(y/2), n);

    %if even square the result
    if (mod(y, 2) == 0)
        result = mod(z*z, n);
        return;
    else
        %odd so square the result & multiply by itself
        result = mod(x*z*z, n);
        return;
    end
end

Found the latter one here: https://gist.github.com/ttezel/4635562