clear all 
clear off
y=zeros(1,150);
for n=1:150        % h= 1/900 之後會出現 round-off
    h=1/(10*n);    % 注意裡面多乘了10，h=1/10,1/20 ...
    p1=log2((N(h)-N(h/2))/(N(h/2)-N(h/4)));
    fprintf('n=%2d  , %1.7f\n',n*10,p1)
    y(n)=p1;
end
loglog(1./[1:150]/10,y)
title('hw3-4-2-12extra-find p1 numerically')
xlabel('h')
ylabel('p1')
grid on

% log10(2)=0.301029995663981,也就是說 p1 約等於 2