Try this:
kx = linspace(-0.5,0.5,15);
ky = linspace(-0.5,0.5,15);
[KX,KY] = meshgrid(kx,kx);
alpha_R = 0.18851;
M = length(kx);
% SPIN MATRICES
sigma_x = [0,1;1,0]; % sigma_x
sigma_y = [0,-1j;1j,0]; %sigma_y
I = [1,0;0,1]; % Identity
%Pre-allocation
E = nan(M,M,2);
Psi = nan(2,1,M,M); % BL's change here
for i = 1: length(kx)
for j = 1: length(ky)
kx_t = kx(i);
ky_t = ky(j);
eps_k = kx_t^2 + ky_t^2; % epsilon_k
%hamiltonian
H = eps_k * I + alpha_R * (sigma_x .* ky_t - sigma_y .* kx_t);
E(i,j,:) = eig(H);
[V,D] = eig(H);
Psi(:,:,i,j) = V(:,1); % BL's change here
end
end
% Psi_prime = pagectranspose(Psi); % BL's change here
%%%NO LONGER ERROR HERE
expx = pagemtimes(pagemtimes(Psi,'ctranspose',sigma_x,'none'),Psi); % BL's change here
expx = reshape(expx,[M M]); % BL's change here
expy = pagemtimes(pagemtimes(Psi,'ctranspose',sigma_y,'none'),Psi); % BL's change here
expy = reshape(expy,[M M]); % BL's change here
% expy
quiver(KX,KY,expx,expy)
