 gp2d(gpmodel, m, s)
%% Code
input = gpmodel.inputs;  target = gpmodel.targets; X = gpmodel.hyp;

if nargout < 4; [M, S, V] = gp2(gpmodel, m, s); return; end

persistent K iK oldX oldIn oldOut beta oldn;           % cache some variables
D = size(input,2);         % number of examples and dimension of input space
[n, E] = size(target);               % number of examples and number of outputs
X = reshape(X, D+2, E);

% 1) If necessary, re-compute cached variables
if numel(X) ~= numel(oldX) || isempty(iK) ||  