This solution uses polyfit and a log transformation, converting the power-law into a linear equation. There might by an easier way. Perhaps someone else will weigh in with an alternate solution.
p = [40 50 60 80 100];
r = [230.88 243.5 268.34 278.92 280];
% transform r = ap^n into log(r) = log(a) + n log(p)
pLog = log(p);
rLog = log(r);
% do the model fitting and get coefficients
pf = polyfit(pLog, rLog, 1);
b = pf(1);
a = exp(pf(2));
% get squared correlation coefficient
CM = corrcoef(pLog, rLog);
R2 = CM(1,2)^2;
% x/y data for power-law curve from x = 0 to x = 150
x = linspace(0,150);
y = a * x.^b;
% plot curve and scatter data
p1 = plot(x, y);
hold on;
s1 = scatter(p, r, 'filled');
ax = gca;
ax.XLim = [0 150];
ax.YLim = [0 350];
ax.XLabel.String = 'p';
ax.YLabel.String = 'r';
% print equation and line to curve
s = sprintf('%s = %5.3f%s^{%.4f}\n%s^2 = %.4f', '{\itr}', ...
a, '{\itp}', b, '{\itR}', R2);
text(50, 150, s);
line([40, 30], [175, a * 30^b], 'color', 'k');
