function[a]=cal(a,b,v)%a,b表示區間，v是精度

i=1;

x=(a+b)/2;

a=[ix];

t=x-(x^3-x-1)/(3*x^2-1);%迭代函式

while(abs(t-x)>v)

i=i+1;

x=t;

a=[a;ix];

t=x-(x^3-x-1)/(3*x^2-1);%迭代函式

end

a=[a;i+1t];

end

執行結果：

>>formatlong;

>>cal(1,2,0.00001)

ans=

1.0000000000000001.500000000000000

2.0000000000000001.347826086956522

3.0000000000000001.325200398950907

4.0000000000000001.324718173999054

5.0000000000000001.324717957244790

主程式：

function [k,x,wuca,yx] = newton(x0,tol)

k=1;

yx1=fun(x0);

yx2=fun1(x0);

x1=x0-yx1/yx2;

while abs(x1-x0)>tol

x0=x1;

yx1=fun(x0);

yx2=fun1(x0);

k=k+1;

x1=x1-yx1/yx2;

end

k;

x=x1;

wuca=abs(x1-x0)/2;

yx=fun(x);

end

分程式1：

function y1=fun(x)

y1=sqrt(x^2+1)-tan(x);

end

分程式2：

function y2=fun1(x)

%函式fun(x)的導數

y2=x/(sqrt(x^2+1))-1/((cos(x))^2);

end

結果：

[k,x,wuca,yx] = newton(-1.2,10^-5)

k =8

x =0.9415

wuca =4.5712e-08

yx =-3.1530e-14

[k,x,wuca,yx] = newton(2.0,10^-5)

k =243

x =NaN

wuca =NaN

yx =NaN