Costas Loop oscillates in ModelSim while locks in MATLAB!

Come on Mr Kazem! The model that we spent days discussing it! We finally said yuppie remember?! Anyway here its is:

% Siraj Muhammad
% 25/3/2015
% BPSK Demodulator
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load RRC.mat
fc = 0.05000001;
phase_offset = pi/7;
N = length(rf_signal);
bb = zeros(1, N);
bb_f = zeros(1, N);
I_f = zeros(1, N);
Q_f = zeros(1, N);
I_r = zeros(1, N);
Q_r = zeros(1, N);
error = zeros(1, N);
PhErr = zeros(1, N);
error_integral = zeros(1, N);
% Loop Filter Coefficients
% BW = 200;  % Hz
% loop_theta = 2*pi*BW;
% C1 = 4*(loop_theta)^2/(1+sqrt(2)*loop_theta+loop_theta^2);
% C2 = 2*sqrt(2)*loop_theta/(1+sqrt(2)*loop_theta+loop_theta^2);
C1 = 2^2;
C2 = 1/2^8;
for i = 2:N
    % Downconverting to Baseband
    bb(i) = rf_signal(i).*exp(j*2*pi*fc*i+j*phase_offset).*exp(-j*PhErr(i-1));
    % Filtering
    bb_f = filter(RRC, bb(1:i));
    I_f = real(bb_f);
    Q_f = imag(bb_f);
    
    % Error
    error(i) = I_f(i).*Q_f(i);
    
    % Loop Filter
    error_integral(i) = error(i).*C1 + error_integral(i-1);
    PhErr(i) = error(i).*C2 + error_integral(i);
end
figure; subplot 321; plot(real(bb)); title('I Channel')
subplot 322; plot(imag(bb)); title('Q Channel')
subplot 323; plot(I_f); title('I Channel Filter')
subplot 324; plot(Q_f); title('Q Channel Filter')
subplot 325; plot(error); title('Phase Error');
subplot 326; plot(PhErr); title('Loop Filter');

OK, that's great... let's stick to modelsim now. so I should expect oscillations in the loop filter output. I believe convergence time may be long and ModelSim is going to consume even longer time to complete simulation.

I will carry on tests tomorrow and let you know... thanks for help Mr Kazem, I am so grateful to you.