numpoints = 1024;

Y1 = fft(noise, numpoints);
Y2 = fft(Vout, numpoints);

Pyy1 = Y1 .* conj(Y1) / numpoints;
Pyy2 = Y2 .* conj(Y2) / numpoints;

f = 10 * 1000 * (0 : floor(numpoints / 2)) / numpoints;
%   1/dt remember ms timescale on some models

loglog(f, Pyy1(1 : floor(numpoints / 2) + 1), f, Pyy2(1 : floor(numpoints / 2) + 1));
xlabel('Frequency (Hz)');
legend('I_{inj}', 'V')
title(strcat('Mean V= ', num2str(round(mean(Vout)))));