Here is the way I would code it. However, because your objective function is just the KL-divergence, it is trivial to see that the minimum occurs at p=pr. In particular, your loop should produce the same result at every iteration.
function [p,fval] = MC_NT(p0,opts,N)
if nargin < 2
opts = optimoptions('fmincon','Algorithm','interior-point');
end
M=length(p0);
p0=p0(:);
pr=0.25*ones(M,1);
p=nan(M,N);
fval=nan(1,N);
for i=1:N
fun=@(p) sum(p.*log(p./pr));
[p(:,i),fval(i)] = fmincon(fun,p0(:),[],[],Aeq1,beq1,[],[],[],opts);
pr=p(:,i);
end
