Volterra Kernel Convolution Method
Yes the MATLAB code is correct. Be careful though, the bandwidth of the signal squared is twice that of the signal itself, which may lead to aliasing if the sampling frequency is too low compared to the signal bandwidth. This can be remedied by properly resampling the signal to a higher sampling frequency with a lowpass filter, then squaring, low-pass filtering, and decimating back to the original rate.

Another pitfall if you are doing everything by yourself is that the impulse-response-like partitions of the deconvolved sweep response are not directly the Volterra kernels. The rest of this answer will discuss how to get usable kernels from the partitions corresponding to the harmonics.

First, you may need to phase shift each partition with a constant phase shift, depending on what kind of sweep you used. You can test if there is any phase shift needed (and test if you do the phase shift correctly), by using as a test system a memoryless non-linearity that is the sum of a number of early Chebyshev polynomials of the first kind, as this will generate for a unit amplitude cosine input a number of harmonics in the same phase (explanation further down in this answer). If I understood the following paper correctly, for certain kinds of sweeps there is no need for phase shifts:

Antonín Novák, Laurent Simon, František Kadlec, Pierrick Lotton. Nonlinear System Identification Using Exponential Swept-Sine Signal. IEEE Transactions on Instrumentation and Measurement Volume: 59, Issue: 8, Aug. 2010.

Then each one-dimensional Volterra kernel can be calculated as the weighted sum of the (phase-shifted) partitions, with the weights obtained from solving a simple linear system. This is covered in Angelo Farina's papers, but at least some of them complicate the derivation by using sines rather than cosines, which gets into the earlier phase shift territory. Cosines are simpler because they have real frequency-domain representations implying no phase shift. Here is an Octave script (using cosines) that solves the weights given the number of harmonics:

N = 3; #Number of harmonics including the fundamental (you can touch this)
M = 2*N+1; #Matrix height, must fit positive, zero and negative frequencies
x = zeros(M, N+1);
# From cosines to complex exponentials. Negative indexing by shifted columns
x(1, 1) = 1; #DC bias
x(3, 2) = 0.5; #This is negative fundamental
x(1, 2) = 0.5; #This is positive fundamental
for k = 3:N+1
  x(:, k) = conv(x(:, k-1), x(:, 2))(1:M); #Time domain exponent k
endfor
for k = 1:N+1
  x(:, k) = shift(x(:, k), 1-k); #Shift back columns
endfor
x = x(1:N+1, :); #Go back to cosines
x(2:N+1,:) .*= 2 #
inv(x) #Solve system


The output is:

x =
   1.00000   0.00000   0.50000   0.00000
   0.00000   1.00000   0.00000   0.75000
   0.00000   0.00000   0.50000   0.00000
   0.00000   0.00000   0.00000   0.25000

ans =
   1  -0  -1  -0
   0   1  -0  -3
   0   0   2  -0
   0   0   0   4


The answer ans is the inverse of matrix x. The answer means that if you have three harmonics including the fundamental, then you go from the partitions $h_n$ corresponding to harmonics $n \in \{1,\, 2,\, 3\}$ to kernels $k_m$ corresponding to exponents $m \in \{1,\, 2,\, 3\}$ by:

$$\left\{\begin{eqnarray}k_0 &amp;=&amp; h_0-h_2\\
k_1 &amp;=&amp; h_1 -3h_3\\
k_2 &amp;=&amp; 2 h_2\\
k_3 &amp;=&amp; 4 h_3\end{eqnarray}\right.$$

This is the solution of:

$$\begin{align}&amp;k_0 + k_1 \cos(t) + k_2 \cos^2(t) + k_3 \cos^3(t) = h_0 + h_1 \cos(t) + h_2 \cos(2t) + h_3 \cos (3t)\\
\Rightarrow &amp; k_0 + \frac{k_2}{2} - h_0 + \left(k_1 + \frac{3k_3}{4} -h_1\right)\cos(t) + \left(\frac{k_2}{2} - h_2\right)\cos(2t) + \left(\frac{k_3}{4} - h_3\right)\cos(3t) = 0,\end{align}$$

where $h_0$ is a zero frequency bias that is not determined by the sweep method and is normally assumed zero, $h_0 = 0,$ and $k_0$ is a term that eliminates the bias. I'm not sure how to use it, probably just add to each output sample the sum of samples of $k_0.$

The solution matrices for smaller N are contained in the solution matrices for larger N. The matrices follow the On-line Encyclopedia of Integer Sequences (OEIS) sequence A053120. OEIS gives simpler ways to calculate the weights, and also hints to a link to Chebyshev polynomials of the first kind. They have coefficients that are identical to the weights. The link is actually such that if you use as the memoryless non-linear elements not powers but Chebyshev polynomials of the first kind, $T_n(x)$, then the kernels to use are directly the (phase-shifted) partials corresponding to the harmonics. This comes from the property $T_n\big(\cos(t)\big)=\cos(nt).$ This is also covered in:

Antonin Novak, Laurent Simon, Pierrick Lotton, Joel Gilbert. Chebyshev Model and Synchronized Swept Sine Method in Nonlinear Audio Effect Modeling. Proc. of the 13th Int. Conference on Digital Audio Effects (DAFx-10), Graz, Austria, September 6-10, 2010.

The N = 7 solution from the script agrees with:

Marc R´ebillat, Romain Hennequin, Etienne Corteel, Brian Katz. Prediction of harmonic distortion generated by electro-dynamic loudspeakers using cascade of Hammerstein models. 128th Convention of the Audio Engineering Society, May 2010, London, United Kingdom. pp.7993, 2010. &lt;hal-00619357&gt;