function [Lambda, Phi] = ComputeKL1D(c, N)

% Mesh
h = 1/N;
x = 0:h:1;

% Matrices of x- and y-coordinates
[Y,X] = meshgrid(x,x);

% Pointwise evaluation of covariance function
tic;
C = c(X,Y);
toc;

% 1D-Trapezoidal rule weight vector
w = h*[0.5, ones(1,N-1), 0.5];
% 2D-Trapezoidal rule weight matrix
e = ones(1, N+1);
W = (w'*e) .* (e'*w);

% Mass matrix
M = h/3 * ( diag([1, 2*ones(1,N-1),1]) + diag(0.5*ones(1,N),1) + diag(0.5*ones(1,N),-1));

% Solve generalized eigenvalue problem
tic;
[V,D] = eig(C.*W, M);
toc;

% Sorting eigenvalues in descending order
d = diag(D);
[d, inds] = sort(d, 'descend');
D = diag(d);
V = V(:,inds);

% Norming eigenfunctions
for i = 1:N+1,
   phi = V(:,i);
   normPhi = sqrt(phi' * M * phi);
   V(:,i) =  V(:,i)/normPhi;
end

Lambda = diag(D);
Phi = V;