Simulating BPSK Costas Loop in Matlab
I figured it out. It turned out that I don't really need that 'Phase Accumulator', and I should reverse the control signal in the NCO. So, I use -PhErr(i-1) as a control signal to NCO. Here is the modified code:

% 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 = 100;  % 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);

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(i) = real(bb_f(i));
    Q_f(i) = imag(bb_f(i));

    % 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');


I added a variable phase_offset to test the loop with different offsets, and I can drift the frequency a little bit and the loop will lock very nicely.

I hope this helps somebody as I couldn't find any MATLAB code for Costas Loop online.