clear 

trueroot=1.1;
Niterations=30;

% root of multiplicity 3
f=inline('2.7951-8.954*x+10.56*x.^2-5.4*x.^3+x.^4');
fp=inline('-8.954+10.56*2*x-5.4*3*x.^2+4*x.^3');

x=0:.01:2;
figure;plot(x,f(x),'r',x,zeros(size(x)),'k') % zeros

% linear convergence -- very slow
x=zeros(Niterations+1,1);
x(1)=1;
for n=1:Niterations
    x(n+1)=x(n)-f(x(n))/fp(x(n));
end
alpha1=abs(x-trueroot);
% alpha(2:end)./alpha(1:end-1).^1

figure;
subplot(1,2,1),plot(log(alpha1),'k')
subplot(1,2,2),plot(alpha1(2:end)./alpha1(1:end-1).^1,'k')

% restore quadratic convergence
x=zeros(Niterations+1,1);
x(1)=1;
for n=1:Niterations
    x(n+1)=x(n)-3*f(x(n))/fp(x(n)); % root of multiplicity 3
end
alpha2=abs(x-trueroot);
% alpha(2:end)./alpha(1:end-1).^1

figure;
subplot(1,2,1),plot(log(alpha2),'r')
subplot(1,2,2),plot(alpha2(2:end)./alpha2(1:end-1).^2,'r') % gets stuck close to root

% better: use 2nd derivative in place of f(x)
f=inline('21.12-32.4*x+12*x.^2');
fp=inline('-32.4+12*2*x');

% x=0:.01:2;
% figure;plot(x,f(x),'r',x,zeros(size(x)),'k') % zeros
x=zeros(Niterations+1,1);
x(1)=1;
for n=1:Niterations
    x(n+1)=x(n)-f(x(n))/fp(x(n));
end
alpha3=abs(x-trueroot);
% alpha(2:end)./alpha(1:end-1).^2 % converges rapidly to correct answer, within machine tolerance

figure;
subplot(1,2,1),hold on; plot(log(alpha3),'b')
subplot(1,2,2),hold on; plot(alpha3(2:end)./alpha3(1:end-1).^2,'b')

figure;
plot(0:Niterations,alpha1,'k',0:Niterations,alpha2,'r',0:Niterations,alpha3,'b')
legend('Newton','NewtonMultRoot','NewtonPrime')