Why MATLAB fft cos makes imaginary parts?
Those $10^{-15}$ terms after the FFT are due to numerical roundoff errors... Apart from that, you should properly frame the cosine into a fully periodic length, as it seems you did, then you can expect the imaginary part to be as close to zero as the numerical format allows...
The following line of OCTAVE / MATLAB code shows the effect of frame length (aka aperture size) on the FFT computation:
L = 16   ; period of cosine
N = 128  ; frame length of sample block (integer multiple of L)
figure,stem(imag(fft(cos(2*pi*[0:N-1]./L),N))); title('imaginary part of the FFT');


Now let us change the period of the cosine so that the frame does not include an exact number of periods:
L = 15  ; period of cosine
N = 128  ; frame length of sample block (non-integer multiple of L)
figure,stem(imag(fft(cos(2*pi*[0:N-1]./L),N))); title('imaginary part of the FFT');


Now if your result display more deviations than those found on the first plot, then your frame length is probably the problem. If you, on the other hand, observe very close to zero (probably randomly distributed) values as in the first plot, then that's probably due to the numerical roundoff errors involved in the FFT computation.