%當f(x)= (10*x.^3+80*x.^2+6*x+a0)*(x-p*)時
[newton2.m]  
function sa = newton2(x0,epsilon,max)
% NEWTON calculate the newton method approximation
% x0 is the initial value 
% epsilon is the tolerance for x0 
% max is the maximum number of iterations
syms x;
syms pstar;
syms q;
syms a0;
q=1000;
pstar=rand(1,1000); % p*
step=zeros(1,1000); %迭代步數
order=zeros(3,1000); %收斂階數
g=inline('-10*x.^3-80*x.^2-6*x');
df=inline('10*x.^3+80*x.^2+6*x+a0+(30*x^2+160*x+6)*(x-pstar)');
U=zeros(1,max);
U(1)=x0;
for i=1:1000
    a0=feval(g,pstar(i));    
    f=inline('(10*x.^3+80*x.^2+6*x+a0)*(x-pstar)');
    for j=1:max
        x1=x0-feval(f,a0,pstar(i),x0)/feval(df,a0,pstar(i),x0);
        U(j+1)=x1;
        err=abs(x1-x0);
        x0=x1;
       if err<epsilon
           break;
       end
    end
    for k=1:j
        Us(k)=log(abs(U(k)-pstar(i)));
    end
    X=Us(1:j-1);
    Y=Us(2:j); 
    if  j==max    %不收斂
        step(i)=0;
        order(:,i)=0;
        q=q-1; 
    else if j==2 %收斂
            X=Us(1:j);
            Y=Us(2:j+1);            
    else          %收斂
        step(i)=j; 
        xmean=mean(X);
        ymean=mean(Y);
        Gs=X-xmean;
        G=zeros(1,j-1);
        for h=1:3
             for d=1:j-1
                 G(d)=Gs(d)*(d^(h-1));
             end
            sumx2=G*(X-xmean)';
            sumxy=(Y-ymean)*G';
            order(h,i)=sumxy/sumx2; 
        end
        end
    end    
end 
q
v=zeros(1,3);
v=sum((order)'); 
averagesteps=sum(step)/q
averageorder1=v(1)/q
averageorder2=v(2)/q
averageorder3=v(3)/q