DouKaku?

function r = doubleperfectnumber(rmax)
% Return every integer, the double of which is equal to the sum of its divisors
% (ja.doukaku.org Q25).
[m n] = meshgrid(uint16(1:rmax-1),uint16(1:rmax));
m(find(m(:)>=n(:))) = 0;
m(find(mod(n,m)~=0)) = 0;
r = find(sum(m,2)==n(:,1)*2);
% plot(1:rmax,sum(m,2),1:rmax,n(:,1));