function pt(m, k)
% test all up to  k prime, all powers from 2 to m.
% fuind factors of p^n - 1, and see whether all are less than p

ps = primes(k);

for p=ps
    disp(p)
    for n = 3:m
        f = factor(p^n - 1);
        b = (f <= p);
        if (all(b))
            'found'
            disp(p)
            disp(n)
            disp(factor(p^n - 1))
        end
    end
end

> pt(7, 500)
     7     4
    41     4
    43     4
    47     4
    67     3
    73     4
    79     3
    83     4
   137     3
   149     3
   157     4
   163     3
   173     4
   181     3
   191     3
   191     4