%Program Newton-Rhapson Iteration
%Computes approximate solution of f(x)=0.
%Input: function f, its derivative df, initial guess x0, and 
%       number of steps k.
%Output: Approximate solution xc and function value at fxc=f(xc).
function [xc,fxc]=newton(f,df,x0,k)
x(1)=x0;

for s=1:k
  x(s+1)=x(s)-f(x(s))/df(x(s));
end
xc=x(k+1); 
fxc=f(xc);
___________________________________________________________________________
% Solving x^3-3x+2 around x=-2 -single root
% Qudratic Convegence 
clear;
format long
iter=6;
f=inline('x^3-3*x+2');
df=inline('3*x^2-3');
A=zeros(iter,3);
for j=1:iter;
[x,y]=newton(f,df,-1.8,j);
A(j,1)=x;
A(j,2)=y;
A(j,3)=abs(-2-A(j,1));
end
B=zeros(iter-1,1);
for j=1:iter-1
    B(j)=A(j+1,3)/(A(j,3).^2);
end
  disp(' ')    %  a blank line
disp(' k           xk                      f(xk)           e(n+1)/en^2 ')
disp('--------------------------------------------------------------------------')
for k=1:iter-1
  disp(sprintf('%2.0f    %7.15e    %7.15e  %7.3e',k,A(k,1),A(k,2),B(k)))
end
 k           xk                      f(xk)           e(n+1)/en^2 
--------------------------------------------------------------------------
 1    -2.033333333333333e+000    -3.067037037037022e-001  6.452e-001
 2    -2.000716845878136e+000    -6.454696479666566e-003  6.662e-001
 3    -2.000000342333275e+000    -3.081000178539739e-006  6.669e-001
 4    -2.000000000000078e+000    -7.034373084024992e-013  0.000e+000
 5    -2.000000000000000e+000    0.000000000000000e+000      NaN
 
%Solving x^3-3x+2 around x=1 -double root
% Linaer convergence 
clear;
format long
iter=15;
f=inline('x^3-3*x+2');
df=inline('3*x^2-3');
A=zeros(iter,3);
for j=1:iter;
[x,y]=newton(f,df,0.8,j);
A(j,1)=x;
A(j,2)=y;
A(j,3)=abs(1-A(j,1));
end
B=zeros(iter-1,1);
for j=1:iter-1
    B(j)=A(j+1,3)/(A(j,3));
end
  disp(' ')    %  a blank line
disp(' k           xk                      f(xk)           e(n+1)/en ')
disp('--------------------------------------------------------------------------')
for k=1:iter-1
  disp(sprintf('%2.0f    %7.15e    %7.15e  %7.3e',k,A(k,1),A(k,2),B(k)))
end
 k           xk                      f(xk)           e(n+1)/en 
--------------------------------------------------------------------------
 1    9.037037037037035e-001    2.692597673118957e-002  4.916e-001
 2    9.526636883316521e-001    6.616111483743081e-003  4.960e-001
 3    9.765230979913695e-001    1.640555138765043e-003  4.980e-001
 4    9.883080249622006e-001    4.085085212048867e-004  4.990e-001
 5    9.941654713260555e-001    1.019265571211570e-004  4.995e-001
 6    9.970855807733605e-001    2.545676367660121e-005  4.998e-001
 7    9.985434992395640e-001    6.361093582984978e-006  4.999e-001
 8    9.992719265315112e-001    1.589886981623323e-006  4.999e-001
 9    9.996360074561469e-001    3.974234905612661e-007  5.000e-001
10    9.998180147710556e-001    9.934984390724821e-008  5.000e-001
11    9.999090101459809e-001    2.483670691333373e-008  5.000e-001
12    9.999545057622717e-001    6.209083025510154e-009  5.000e-001
13    9.999772530543009e-001    1.552258988013477e-009  5.000e-001
14    9.999886265713093e-001    3.880629151353787e-010  5.000e-001
%__________________________________________________________________________
%Modified version
%Modified Newton-Rhapson Iteration
%Computes approximate solution of f(x)=0 at double root.
%Input: inline function f, its derivative df, initial guess x0, and 
%       number of steps k
%Output: Approximate solution xc and function value and f(xc).
function [xc,fxc]=newton1(f,df,x0,k)
x(1)=x0;
for s=1:k
  x(s+1)=x(s)-(2*f(x(s)))/df(x(s));
end
xc=x(k+1); 
fxc=f(xc);
___________________________________________________________________________
 %Solving x^3-3x+2=0 around x=1 -double root
% Quadratic convergence 
clear;
format long
iter=6;
f=inline('x^3-3*x+2');
df=inline('3*x^2-3');
A=zeros(iter,3);
for j=1:iter;
[x,y]=newton1(f,df,0.8,j);
A(j,1)=x;
A(j,2)=y;
A(j,3)=abs(-2-A(j,1));
end
B=zeros(iter-1,1);
for j=1:iter-1
    B(j)=A(j+1,3)/(A(j,3).^2);
end
  disp(' ')    %  a blank line
disp(' k           xk                      f(xk)           e(n+1)/en^2 ')
disp('--------------------------------------------------------------------------')
for k=1:iter-1
  disp(sprintf('%2.0f    %7.15e    %7.15e  %7.3e',k,A(k,1),A(k,2),B(k)))
end
 
 k           xk                      f(xk)           e(n+1)/en^2 
--------------------------------------------------------------------------
 1    1.007407407407407e+000    1.650154956054362e-004  3.317e-001
 2    1.000009111202218e+000    2.490425643486560e-010  3.333e-001
 3    1.000000000021503e+000    0.000000000000000e+000  3.333e-001
 4    1.000000000021503e+000    0.000000000000000e+000  3.333e-001
 5    1.000000000021503e+000    0.000000000000000e+000  3.333e-001