You are quite right. BPSK signal, sampled at Nyquist rate, is indeed just a sequence of independently generated -1 and 1 - so you got that right. For convenience, let us denote this discrete time sequence by $x_n$ ($n$ denotes time domain).
Your second observation was that the frequency transform was taking up the entire spectrum, and had peaks at random frequency locations. That is expected as well. Below, I explain why this happens and how to estimate the power spectrum.
Consider the time sequence $x_n$. The power spectrum gives you how much energy it occupies in different frequencies. To calculate the power spectrum, we start with the time domain autocorrelation and take its Fourier transform.
Since the sequence $x_n$ is independent, its autocorrelation $R(\tau)$ has energy only at the zeroth term i.e.
$$R(\tau) = E[x_n x_{n+\tau}^*] = \left\{ \begin{array}{cc} 1 &if \tau=0 \\ 0 & if \tau \neq 0 \end{array} \right.$$
Denoting the power spectrum of $x_n$ by $S(f)$, we have the following
$$ S(f) = \sum_{k=-\infty}^{+\infty} R(k) e^{-j2\pi f k T_s} \\
= R(0) e^{-j2\pi f \times 0 \times T_s} \\
= 1 $$
Here, $T_s$ is the sampling period associated with the discrete time sequence $x_n$, and the frequency range is $f \in [-1/{2T_s}, 1/{2T_s}]$. The second step in the equation above makes use of $R(k)=0$ if $k\neq 0$.
So, we expect $x_n$ to have a flat frequency response; i.e. all frequencies are expected to have equal energy. We do not expect energy to peak at DC for this input sequence.
However, in your Matlab code, you observed randomly located peaks in the frequency transform. This happens because you took the FFT of $\{x_1,\cdots,x_{1000}\}$; i.e. this is because you are looking at the frequency response of a specific instance of a random signal.
The power spectrum does not really have these peaks. To get the power spectrum, you need to do some additional steps. What you need is the average energy of the frequency response. I suggest adding an outer loop to your code, as below, and averaging the energies you measure in each frequency.
num_iter = 10000;
sample_len_per_iter = 1000;
power_spectrum = zeros(1,sample_len_per_iter);
for iter_index = 1:num_iter
bits=randi([0 1],1,sample_len_per_iter);
bpsk_mod=2*bits-1;
fft_bpsk=fft(bpsk_mod);
power_spectrum = power_spectrum + abs(fft_bpsk.^2); % accumulate power
end
power_spectrum = power_spectrum / num_iter; % average out
f=[-sample_len_per_iter/2:sample_len_per_iter/2-1];
plot(f,power_spectrum)
The above should give you a much more flat response.
In case you are curious, the fft output of the specific instance of BPSK sequence is (or converges to) an i.i.d. complex Gaussian sequence. This is a consequence of the central limit theorem. This is why you observed randomly located peaks.
In a practical system, this signal is up-sampled and filtered by a baseband pulse shape so that the bandwidth actually occupied meets some pre-specified requirement of RF emissions.
Edit about the power scaling:
The fft() function in Matlab includes a scaling of $1/\sqrt{L}$ where L is the fft length i.e. Length of fft_bpsk. This, after taking magnitude square, gives the desired power scaling of $1/L$.
$$X_k = \frac 1 {\sqrt L} \sum_{n=1}^{L} e^{-j2\pi \frac{(n-1)(k-1)}{L}} x_n\\ E[|X_k|^2] = \frac 1 L \sum_{n=1}^{L} \sum_{m=1}^{L} e^{-j2\pi \frac{(n-m)(k-1)}{L}} E[x_n x_m^*] \\ = \frac 1 L L R(0) = 1 $$