How can I compute a time-frequency cross-spectrum in MATLAB?
I don't know of a function that does all of what you ask, but it is easy enough to write.  They key step is to create short-time integrations on which the FFTs may be performed.  As you say, cpsd then averages all these STIs, but you can write your own to skip that step.

The key function to do this is $y = buffer (x,fftlen,overlap)$.  You input your FFT processing parameters, similar to cpsd, and get a matrix of overlapping time series.  I apply a Hanning window to each STI before performing FFT. Then you can simply apply $S = FFT(y)$ to get a matrix of overlapped spectra.

So it would look something like this:

% Assuming x is longer than fftlen:
overlap = fftlen/2; % 50% overlap
win = hanning(fftlen);
X = buffer(x,fftlen,fftlen/2,'nodelay');  % Matrix of overlapping STIs
numSTIs = size(X,2);
winX = X.*win(:,ones(1,numSTIs); % Time-domain windowed STIs
S = fft(winX,fftlen,1)/fftlen;  % Double-Sided Complex Spectrum Matrix
SdB = 20*log10(2*abs(S(1:fftlen/2+1,:)));  % Log Scale Single-Sided Real Spectrum Matrix


You can then calculate your auto-spectrum using either of these matices.  Of course you would need an x2 to repeat and calculate a cross-spectrum matrix.

for i = 1:numSTIs
    Sx1x2(:,:,i) = S1(:,i)'*S2(:,i);
end


(If anyone knows of a way to do this without a for loop, please let me know!)