I am new to Matlab and currently trying to pass an online course on it. The task is to compute Feigenbauim's constant, delta, which is the limiting ratio of each bifurcation interval to the next between every period doubling, of a one-parameter logistic map, to an arbitrary degree of accuracy. Here is my script:

% Compute the Feigenbaum delta
% Store approximate values in the row vector delta for assessment, where length(delta)= num_doublings and 
% delta(2:num_doublings) are computed from the algorithm described in Lectures 21-23.
num_doublings=11; delta=zeros(1,num_doublings); delta(1)=5;
% Write your code here

m(1)=2;
m(2)=1+sqrt(5);


for i=2:num_doublings
    x0 = 1/2; 
    dx0 = 0;
    period =2^(i-1);
    m(i+1)=m(i)+(m(i)-m(i-1))/delta(i-1);
   
    for j=1:100
        for k=1:period
            x = m(i+1)*x0*(1-x0);
            x_p = x0*(1-x0) + m(i+1)*dx0*(1-2*x0);
            
            x0 = x; dx0 = x_p;
        end
        u(j)= m(i+1)-(x0-0.5)/dx0;
    end
    m(i+1)=u(j);
    delta(i)=(m(i)-m(i-1))/(m(i+1)-m(i));
end

% Output your results
fprintf('n        delta(n)\n');
for n=1:num_doublings
    fprintf('%2g %18.15f\n',n,delta(n));
end

It computes the constant only to the fifth decimal place, and don't know, why if anyone would be willing to give any advice on how to improve it, it would be appreaciated