FFT

Hello,

I have a system of differential equations and I apply a sinusoidal perturbation in current to get the output voltage. Finally, I want to calculate the impedance values and draw the nyquist plot. I started with a simple example of RC circuit and solve the system of equation using ODE45 ( instead of considering the impedance values 1/jwc & R: I started with this simple example to find out how I can buid up the nyquist plot from time-domain signals. my real system is a system of 10 ODEs , don't have the equivalent circuit and want to use the governing equations)

function dx=RCsystem(t,x)
dx=zeros(1,1);
      R=50;      C=6e-5;% f=1;
frequ=0.1;
current=1*sin(2*pi*frequ*t); % i time domai
dx(1)=current/C-x(1)/(R*C); % apply perturbation in current
end 
[t x]=ode45(@RCsystem,tspan,0);

now I have the discrete signal x in time and I want to find the related impedance. If I change the values of frequency and find the corresponding impedance at each frequency ( 1HZ to 1e5 HZ), I can draw the nyquist plot. My question is when I have the time domain signal (discrete), resulted from the integration of differential equations, at different range of frequencies, how can I convert the signals to frequency domain and find the corresponding points in the nyquist plot?

As an example, I tried this with frequ=0.1 in RCsystem function. But, the impedance value is not correct. I appreciate if you help me and tell me where I am doing wrong.

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 1000;                     % Length of signal
t = (0:L-1)*T;                % Time vector
tspan=t;
[t x]=ode45(@RCsystem,tspan,0);
y = x;     % signal from integration
plot(Fs*t(1:50),y(1:50))
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
% Plot single-sided amplitude spectrum.
figure(1)
plot(f,2*abs(Y(1:NFFT/2+1))) 
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
figure(2)
plot(real(Y),imag(Y),'ok')
% find the maximum peak in periodogram and find the corresponding index (I)
% read the value (real and imaginary part) this is one point on nyquist
% plot realated to the specific frequency
maxima=2*abs(Y(1:NFFT/2+1));
[c,I]=max(maxima);
figure(3);plot(real(Y(I)),-imag(Y(I)),'ok')

Thanks Mona