Swap your definitions of m and n (or equivalently use rand(n,m)). Like you said, alpha is 1st dimension, beta is 2nd.
alpha = 1:10;
beta = 1:8;
gamma = 1:4;
lambda = 1:4;
m = length(alpha);
n = length(beta);
Data(:,:,1,1) = rand(m,n); % alpha in first dimension, beta is 2nd dimension, gamma in 3rd and lambda in 4th
Data(:,:,2,1) = rand(m,n);
Data(:,:,3,1) = rand(m,n);
Data(:,:,4,1) = rand(m,n);
Data(:,:,1,2) = rand(m,n);
Data(:,:,2,2) = rand(m,n);
Data(:,:,3,2) = rand(m,n);
Data(:,:,4,2) = rand(m,n);
Data(:,:,1,3) = rand(m,n);
Data(:,:,2,3) = rand(m,n);
Data(:,:,3,3) = rand(m,n);
Data(:,:,4,3) = rand(m,n);
Data(:,:,1,4) = rand(m,n);
Data(:,:,2,4) = rand(m,n);
Data(:,:,3,4) = rand(m,n);
Data(:,:,4,4) = rand(m,n);
idat = interpn(alpha,beta,gamma,lambda,Data,1,2,3,4);
disp(idat)
Edit to show that the interpolated value matches the input table:
Data(1,2,3,4)
