Help on calculating the model

From the estimated loadings it is possible to estimate the model of the data as

Xk = ADkBT,

where Dk is a diagonal matrix holding the kth row of C in its diagonal and Xk is the k’th frontal slab of the three-way array (see figure below). This notation ADkBis quite common in multi-way analysis. It is also quite difficult to see why that is a PARAFAC model!

However, assume that the first mode in your array is the excitation mode and the second is the emission mode. Then A and B will hold the pure excitation and emission spectra.
The third mode (corresponding to different samples) loading or score matrix C will then hold the concentrations. The k‘th row in C will hold the concentrations of the k‘th sample.

How can Xk be modeled? 
Xk is holding the spectra measured in the k‘th sample so that must be the sum of the contribution from the first analyte plus the contribution from the second analyte etc.
The contribution from the first analyte is a1b1T times the concentration of the analyte which is ck1, so the contribution is a1ck1b1T. 

Similar holds for the other analytes, so Xmust be the weighted sum of all the bilinear parts a1b1Ta2b2T, etc. where each bilinear part (factor) is weighted by the concentration. This is exactly what is stated in the expression ADkBT, because the diagonal matrix Dk holds the concentrations.

Try to Think about it, make drawings, and see if it makes sense.

This is equivalent to
X(I×JK) = X = A(C×B)T.

Where x is the so-called Kathri-Rao columnwise Kronecker product.
This product enables one to express the model of the whole data set
(rather than just an individual slab) in terms ordinary algebraic terms.

You may calculate the model of X called directly as
M = A*kr(C,B)';

But then you will get the matricized version of the model which will be of size 5×12261
(in case of using the full data set).

You can convert to the three-way array by 
M = reshape(M,[5 201 61]);

Alternatively you can simply use the function nmodel:
M = nmodel(Factors);