% define the kernel
%k = @(x,y) 1*x'*y;
k= @(x,y) exp(-0.5*(1/0.25^2)*(x-y)'*(x-y));

% sample points
x= (0:.005:1);
n = length(x);

% construct covariance matrix
C = zeros(n,n);
for i=1:n
    for j=1:n
        C(i,j)= k(x(i),x(j));
    end
end

% sample from GP at these points
u = randn(n,1);
[A, S, B]= svd(C);
z = A*sqrt(S)*u;  % z = A S^0.5 u

figure(2); hold on;
plot(x,z,'.-')
axis([0, 1, -2, 2])

