Matlab: Help in understanding if the example for Maximum Likelihood Equalizer is properly functioning
I had a look at your code, and I know why you get 0 BER.

The thing is that you do not add any noise, you should use Eb/N0 to calculate the noise variance, and add a complex noise on your channel output.

The distortion added by "chanOutput = filter(chCoeffs,1,modSignal);" merely adds linear ISI, by convolving the channel impulse response with the input signal. This linear distortion is exactly known to the equalizer since you pass it the exact channel coefficients. This is never the case in real applications.

Even assuming you do know the channel coefficients, you still need to add noise immediately after the channel distortion and before the equalizer, to get a meaningful result of BER over Eb/N0. 

Noise can be generated with randn having unit variance, then you need to multiply this noise by a constant.

I did not spend too much time understanding what you want to verify, but i did correct the script to include noise.

clear
N = 10^4; % number of symbols
N2= 10;   %repeat with same Eb/N0

% Eb_N0_dB = [0:25]; % multiple Eb/N0 values
Eb_N0_dB = (0:10); % multiple Eb/N0 values
nTx = 1;  %SISO 
nRx = 1;

M = 4; %QPSK

for ii = 1:length(Eb_N0_dB)

    %% ----- ADDED
    EbN0 = Eb_N0_dB(ii);
    SNR_linear  = 10^(EbN0/10)*log2(M);
    SigPWR      = 1; %unit circle symbols
    NoisePWR    = SigPWR/SNR_linear;
    NoiseFactor = sqrt(NoisePWR/2); %half power complex half imaginary

    hMod = comm.QPSKModulator(0,'SymbolMapping','Binary');
    hDemod = comm.QPSKDemodulator(0,'SymbolMapping','Binary');
 % Channel coefficients
    chCoeffs = [.986; .845; .237; .12345+.31i]; 
    hMLSEE = comm.MLSEEqualizer('TracebackDepth',10,...
                      'Channel',chCoeffs, 'Constellation',[1 1i -1 -1i]);
 % Create an error rate calculator
    hError = comm.ErrorRate;
    for n = 1:N2
      data= randi([0 3],N,1);
      modSignal = step(hMod, data);

 % Introduce channel distortion.
      chanOutput = filter(chCoeffs,1,modSignal); 

      %% ----- ADD NOISE !!      
      noiseVec   = NoiseFactor*randn(N,1)+1j*NoiseFactor*randn(N,1);
      chanOutput = chanOutput+noiseVec;


 % Equalize the channel output and demodulate
      eqSignal = step(hMLSEE,chanOutput);
      demodData = step(hDemod,eqSignal);
 % Compute BER
      errorStats = step(hError, data, demodData);
%       nErr=errorStats(1) ;
    end
    fprintf('Error rate = %f\nNumber of errors = %d\n', ...
      errorStats(1), errorStats(2))

    %assume 1 symbol error = 1 bit error (needs Gray coding in real system for this to be true)
    SER(ii) = errorStats(1); % SER ~= BER with undelrying Gray coding
end

%% changed, this is already divided by N*N2
simBer = SER; % simulated ber


EbN0Lin = 10.^(Eb_N0_dB/10);
theoryBer_nRx1 = 0.5.*(1-1*(1+1./EbN0Lin).^(-0.5)); 
p = 1/2 - 1/2*(1+1./EbN0Lin).^(-1/2);
theoryBerMRC_nRx2 = p.^2.*(1+2*(1-p)); 

close all
figure
semilogy(Eb_N0_dB,theoryBer_nRx1,'bp-','LineWidth',2);
hold on
semilogy(Eb_N0_dB,theoryBerMRC_nRx2,'kd-','LineWidth',2);
semilogy(Eb_N0_dB,simBer,'mo-','LineWidth',2);
axis([0 25 10^-5 0.5])
grid on
legend('theory (nTx=1,nRx=1)', 'theory (nTx=1,nRx=2, MRC)', 'sim (nTx=2, nRx=2, ML)');
xlabel('Average Eb/No,dB');
ylabel('Bit Error Rate');
title('BER for BPSK modulation with ML equalizer');