[X1,X2] = meshgrid(-pi:pi/16:+pi, -pi:pi/16:+pi);

Y = sin(X1).*sin(X2) + 0.1*randn(size(X1));

imagesc(Y); drawnow;

x = [X1(:) X2(:)]; y = Y(:);

covfunc = @covSEiso; likfunc = @likGauss; sn = 0.1; hyp.lik = log(sn);

hyp2.cov = [0 ; 0];
hyp2.lik = log(0.1); 
hyp2 = minimize(hyp2, @gp, -100, @infExact, [], covfunc, likfunc, x, y); 
exp(hyp2.lik)
nlml2 = gp(hyp2, @infExact, [], covfunc, likfunc, x, y)

[m s2] = gp(hyp2, @infExact, [], covfunc, likfunc, x, y, x);

m = reshape(m, size(Y));

figure(2); imagesc(m);