function R = correlationfunc(l,x1,x2,type)
if nargin == 3 %'gaussian'
    R = exp(-0.5*pdist2(x1,x2).^2/(l^2));
else
    switch type
        case 'matern52'
            mdist = pdist2(x1,x2)/l;
            m5f = prod(1+sqrt(5).*mdist+(5/3).*mdist.^2,3);
            m5s = exp(-sqrt(5)*sum(mdist,3));
            R = m5f.*m5s;
        case 'gauss'
            R = exp(-0.5*pdist2(x1,x2).^2/(l^2));
    end
end
R(R<eps)=0;
end