I am solving a system of discretized ODEs using $N=256$ grid points and $L = 256$ is the length of the domain. The spatial grid is equally spaced -L/2:1:L/2-1 Using Matlab’s ode23 solver, the resulting solution vector u is a 2177 x 256 matrix in which, per the documentation, each row in u corresponds to a time in the column vector t. That is, u has $256$ components and each row is simply a time slice.

I want to compute the Power Spectral Density $|hat{u}(t)|^2$ plot it against the discretized wavenumbers $kappa$ from $-pi$ to $pi$.

According to the fft(X,n) documentation,

If X is a matrix, then each column is treated as in the vector case.
If X is a vector, then fft(X) returns the Fourier transform of the
vector.

What confuses me is that I have more than one component in my solution vector which differs than the examples I have read. To use the fft command as desired, do I:

1. set n = length(t) = 2177 and use uhat=fft(u,n)? The PSD would then be PSD = uhat.*conj(uhat)/n where we divide by n to normalize.
2. How would I plot the PSD against wavenumbers?