I'm not familiar with the term energy spectral density, you're probably referring to power spectral density. In your case, it is probably sufficient to estimate the PSD using the DFT. Your chirp in the frequency domain is Y = fft. Then to get the PSD, P = abs(Y.^2). Now, you've got to be careful. When we take the DFT using MATLAB, it uses the formula:
Y(k) = sum (from j = 1 to n) X(j)Wn^(j-1)(k-1) (where Wn = e^(-2*pi*i/n) )
(all according to the documentation page)
the frequency of the complex exponential is 'proportional' to k. In reality, the frequency that we know reaches a maximum when k = n/2 + 1 because that is the highest frequency discrete signal (The signal we are mixing is exp(2*pi*i*n/2*j/n) = exp(pi*i*j) = (-1)^j. In visualization, we are used to putting the highest frequency toward the right and the lowest at the middle, but in this case we have the lowest at the left and the highest in the middle. To fix this, use P_shifted = fftshift(P) and then display that.