# 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.