Problem implementing KL expansion with square-exponential kernel: output just looks like Gaussian I am trying to implement the Karhunen Loeve expansion for a 1-D Gaussian random field with a square-exponential kernel. 

Specifically, I know that a Gaussian process has a KL expansion $\hat{U}=\sum_{i-1}^{N_{KL}} \sqrt{\lambda_i}u_i(\mathbf{x})\xi_i(\theta) $, where $\xi_i$ are independent, center-normalized Gaussian variables, and $\lambda_i$ and $u_i$ are the eigenvalue/eigenfunction solutions to the Fredholm equation $\int_\Theta K(\mathbf{x},\mathbf{x}')u_i(\mathbf{x}')d\mathbf{x}'=\lambda_iu_i(\mathbf{x})$.

For the case of a square-exponential kernel of the form $K(\mathbf{x},\mathbf{x}')=\sigma^2\exp(\frac{-(x-x')^2}{2l^2})$, This source gives the solution $\lambda_k=\sqrt{\frac{2a}{A}}B^k$ and $u_k(x)=\exp(-(c-a)x^2)H_k(\sqrt{2c}x)$, where $H_k$ are $k$th order Hermite polynomials, $a^{-1}=4\sigma^2$, $b^{-1}=2l^2$, $c=\sqrt{a^2+2ab}$, $A=a+b+c$, and $B=b/A$. 

I am now trying to apply this to generate a Gaussian process for a domain length of 1000 units, with a grid spacing of 1 unit, and correlation length of 20 units. Because random variables in a KL expansion are drawn from a unit-normal distribution, I interpreted that as $sigma^2=1$. My MATLAB code for this is below:



% define domain:
x = linspace(0,999,1000);

% define statistical properties:
correlation_length = 20;
variance = 1;
mean = 0;

% helper variables for eigenvalues/functions:
a = (4*variance.^2).^-1;
b = (2*correlation_length.^2).^-1;
c = sqrt(a.^2 + 2*a*b);
A = a + b + c;
B = b/A;

% set arbitrary truncation point to 20 eigenfunctions:
Nkl = 20;
% plot 4 realizations: 
U = zeros(4,1000);
for i = 1:4
   for k = 1:Nkl 
       U(i,:) = U(i,:) + sqrt(kth_eval(k, a, A, B)) .* ...
           kth_efunc(k, a, c, x) * normrnd(0,1);
   end
end
figure;
plot(x, U);

% eigen value/vector calculations:
function e_val = kth_eval(k, a, A, B)
    % return scalar value of kth eigenvalue
    e_val = sqrt(2*a/A) * B.^k;
end

function e_func = kth_efunc(k, a, c, x)
    % return 1-D vector of values taken by the kth eigenfunction 
    % over the domain x
    e_func = exp(-(c-a)*x.^2).*hermiteH(k,sqrt(2*c)*x);
end


This gives me an output with the following appearance:



Which seems incorrect to me, as I want to get a process that is statistically similar to the following:



With this code, I can also plot the $k=1,2,3$ eigenfunctions, giving the plot: 



1 has a plot of the first $k=0,1,2$ eigenfunctions, as shown below for $a=1, b=3$:



When I plot the first 3 eigenfunctions with my code with the same a and b, I do get a visually similar results:



However, even this still generates problems for me, as when I attempt to realize a Gaussian process as with the code above and normalized random variables, I seem to just get the same symmetric shape re-scaled by a random factor, rather than random shapes that are statistically similar:



If anyone knows what's causing the discrepancy between the results of my code and what I am expecting, I would be grateful to know.
I've had issue with this -- please see the errata of the reference. The eigenfunctions are not scaled in the original version of the book. This errata should solve it.