%f(x)= 10*x.^3+80*x.^2+6*x+a0時(shí)，
[secant1.m]
function sa=secant(x0,x1,epsilon,max)
% NEWTON calculate the newton method approximation
% x0,x1 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;
m=1000;
pstar=rand(1,m);  %p*
step=zeros(1,m); %迭代步數(shù)
order=zeros(3,m); %收斂階數(shù)
g=inline('-10*x.^3-80*x.^2-6*x');
for i=1:m
    a0=feval(g,pstar(i));    
    f=inline('10*x.^3+80*x.^2+6*x+a0');
    U=zeros(1,max);
    U(1)=x0;
    U(2)=x1;
    for j=1:max
        x2=x1-feval(f,a0,x1)*(x1-x0)/(feval(f,a0,x1)-feval(f,a0,x0));
        U(j+2)=x2;
        err=abs(x2-x1);
        x0=x1;
        x1=x2;
       if err<epsilon
           break;
       end
    end
    for k=1:j+1
        Us(k)=log(abs(U(k)-pstar(i)));
    end
    X=Us(1:j);
    Y=Us(2:j+1);
    if  j==max    %不收斂
        step(i)=0;
        order(:,i)=0;
        q=q-1; 
    else          %不收斂
        step(i)=j;
        xmean=mean(X);
        ymean=mean(Y);
        Gs=X-xmean;
        G=zeros(1,j);
        for h=1:3
            for d=1:j
                G(d)=Gs(d)*(d^(h-1));
            end
            sumx2=G*(X-xmean)';
            sumxy=(Y-ymean)*G';
            order(h,i)=sumxy/sumx2; 
        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