Constellation diagram for QPSK modulation
I think you may be missing some fundamental understanding of QPSK modulation, which was likely the point of the homework exercise so let me hopefully fill in some blanks, from which you can figure out how to code up in your Matlab.

This diagram of a QPSK modulator may be helpful, taken from http://www.masoodkh.com/projects/communications/18-qpsk-modulatordemodulator: 



Note the similarity to your implementation, where the data is converted to +/-1, and multiplied by the sine and cosine of the carrier.  The process of doing this takes your waveform that would otherwise be centered at frequency = 0 (DC) which is the baseband waveform and move it to a carrier frequency (where it is easier to move the signal over great distances). With such a linear modulator, everything about that signal at DC is the same at the carrier in terms of the signal's amplitude and phase, we just change the frequency the signal rides on. So we start with a "DC" carrier and move it to $f_c$ which is what the block diagram is doing and what you are doing in your code. The constellation is the baseband amplitude and phase assignments that are made at DC prior to shifting the signal up to a higher frequency carrier.



We call the baseband data IQ data as it is typically a complex data stream with real and imaginary components, "I" is in-phase and represents the real portion and "Q" is the quadrature-phase and represents the imaginary components. This is certainly the case with QPSK, we take every other bit and assign it to the real axis and the others to the imaginary axis to get I+jQ at baseband.

Your original constellation and how to plot it without the noise should now be very straight forward (plot I+jQ before it gets multiplied by the sine and cosine functions).

To plot the constellation after adding noise seems to me would be very challenging to ask for a homework question, which makes me think there is an easier way than what I proceed to describe below (and hopefully someone else can provide that), but here are my thoughts:

It is not clear if the constellation with noise should be plotted for all samples over the symbol duration, or if you are to get an average of the noise to plot the one sample that a demodulator would determine before making a decision (it would be the latter if the assignment was clear that there should be only one sample per symbol), but let me describe both cases:

Note that the "carrier" is multiplying your baseband signal by $e^{j2\pi f_c t}$ since

$$e^{j2\pi f_c t}= \cos(2\pi f_c t)+ j\sin(2\pi f_c t)$$.

We assign a $j$ to the Q input of the modulator and multiplying by this identity according to Euler matches the block diagram above. The implementation is specifically equal to taking the real part as follows:

$$Re[(I+jQ)(e^{-j2\pi f_c t})] = I\cos(2\pi f_c t)+jQ\sin(2\pi f_c t)$$

Importantly based on the 4 possible dibits on the input ({1,1},{1,-1},{-1,1},{-1,-1}) we assign 4 phases (45°, 135°, 225°,315°) to the carrier. Note that the math and seeing the result is actually much easier using the $e$ instead of sines and cosines, understanding that $Ke^{j\theta}$ is a vector on the complex plane with magnitude $K$ and angle $\theta$.  

Note that multiplying a signal by $e^{j\omega_c t}$ (where $\omega_c=2\pi f_c$) simply shifts the signal in frequency in one direction, depending on the sign of $\omega_c$:

A modulated signal that is varying in amplitude and phase (and in your case of unfiltered QPSK, varying predominantly in phase only) will have the same amplitude and phase variation at baseband (DC) as at any carrier, as represented in the IQ diagrams below. It's just in the case of a carrier, the amplitude and phase is plotted with reference to the carrier. This plot shows how multiplying a signal at the carrier frequency by $e^{-j2 \pi F_o t}$ (where $F_o$ was used in the diagram but is the same as your $f_c$) will move the signal from the carrier back to baseband (DC). 



Ultimately to see your constellation with the noise, you want to move your signal that is at the carrier back to DC. A further complication is that your signal as modulated is completely real (note that I took the real part after multiplying by $e^{j2 \pi f_c}$ when we did the upconversion to make it match your implementation. A complex signal can be single sided (positive frequencies only for example) as we described would happen after multiplying by $e^j2\pi f_c t$ (and note that I specifically multiplied by $e^-j2\pi f_c t$ in order to match your implementation, which means I shifted the frequency to a negative frequency carrier!). Once we take the real part of a complex frequency, what was once a single sided spectrum (if it was, doesn't have to be but will always be the case is the positive spectrum is different from the negative spectrum), will become a double sided spectrum that is complex conjugate symmetric: For all real signals the spectrum will always be complex conjugate symmetric- meaning the negative frequency spectrum will equal the positive frequency spectrum in magnitude and opposite phase.  

Since your signal is real, there will be two complex conjugate symmetric spectrums, one at $+f_c$ and one at $-f_c$. (When I took the real part of what was originally just a spectrum at $-f_c$, the companion spectrum at $+f_c$ is immediately created by the operation of taking the real part).

That said, to move your signal back to DC, simply multiply by $e^{-j2\pi f_c t}$ which will shift both spectrums to the left, and you will get one at DC and another that is at $-2 f_c$. It is possible here that the spectrum will be inverted, as I did not follow the math carefully in my processing described, in that case it would be better to multiply by $e^{-j2\pi f_c t}$, or simply swap I and Q and change the sign of Q which will also invert the spectrum. You then need to filter out the higher frequency signal that will be at twice the carrier frequency in order to get your baseband constellation only. This is the part I am surprised you would have to do as a homework problem; and if you had access to the noise samples at baseband this would be much easier but given the way the noise was added it doesn't appear that you do?

If you needed to plot the constellation with noise as 1 sample per symbol (IQ dibit), then you would simply take the average over a symbol period of the signal after converting back to baseband. Actually, by taking the average, the required low pass filtering will also be done (by nature of averaging is a low pass filter). So such a result for the constellation with noise is actually easier achieved when plotting one sample per symbol. Note too that taking such an average, if synchronized with the symbol period, represents the ideal matched filter for your example of QPSK with a rectangular pulse shape. This is often depicted on diagrams of receiver implementations as an integrate and dump filter (moving average where you accumulate or integrate over the symbol period, dump the result and reset the moving average).