Each grad is a 2x2 matrix; it's not clear from the question what size you want the final grad array to be, so here's one way to put them all in a 2x2xN array:
N = 4; % number of experiments
result = zeros(2*N,1); % pre-allocate result arrays
grad_all = zeros(2,2,N);
dt=0.01;
T=1;
k=[5;1];
y0_all = [2 1 1 4; 1 1 0 1];
for ii = 1:N
y0 = y0_all(:,ii);
[SP,grad]=enzyme(y0,k,dt,T);
result((ii-1)*2+[1 2]) = y0-SP;
grad_all(:,:,ii) = grad;
end
disp(grad_all)
