echo off;
% NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
% To accompany the text:
% NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
% Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
% This free software is complements of the author.

% Algorithm 2.8 (Muller's Method).
% Section 2.5, Aitken's Process & Steffensen's & Muller's Methods, Page 97
echo on; clc; format long; hold off; clear
% This program implements Muller's method.

% Define and store the function f(x) in the M-file f.m
% function y = f(x)
% y = x.^3 - 3.*x + 2;

delete f.m
diary f.m; disp('function y = f(x)');...
disp('y = x.^3 - 3.*x + 2;');...
diary off;

% Remark. f.m and muller.m are used for Algorithm 2.8
f(0); % Test for file f.m
pause % Press any key to see the graph y = f(x).

clc;
% Plot f(x) over the interval [a,b].

a = -2.5;
b = 2.5;
c = -5;
d = 5;
h = (b-a)/100;
X = a:h:b;
Y = f(X);
axis([a b c d]);...
plot(X,Y,'-g');...
hold on;...
plot([a b],[0 0],'b',[0 0],[c d],'b');...
xlabel('x');...
ylabel('y');...
title('Graph of y = f(x).');...
grid;...
axis;...
hold off;...
shg; pause % Press any key to perform Muller's iteration.

clc;
% Place the starting values in p0 p1 and p2

% Place the abscissa tolerance in delta

% Place the ordinate tolerance in epsilon

% Place the number of iterations in max1

p0 = -2.6;
p1 = -2.5;
p2 = -2.4;
delta = 1e-12;
epsilon = 1e-12;
max1 = 10;

[p2,y2,err,P] = muller('f',p0,p1,p2,delta,epsilon,max1);

pause % Press any key for the list of iterations.

clc; clg;
a = -2.65;
b = -1.95;
c = -9;
d = 1;
h = (b-a)/100;
X = a:h:b;
Y = f(X);
max1 = length(P);
n0 = min(7,max1);
X0 = P(1:n0);
Z0 = zeros(1,n0);
axis([a b c d]);...
plot(X,Y,'-g',X0,Z0,'or');...
hold on;...
plot([a b],[0 0],'b',[0 0],[c d],'b');...
xlabel('x');...
ylabel('y');...
title('Graphical analysis for Muller`s method.');...
grid;...
axis;...
hold off;...
shg; pause % Press any key to continue.

J = 1:max1;
Yp = f(P);
points = [J;P;Yp];
Mx1 = 'Iterations for Muller`s method.';
Mx2 = [' k p(k) f(p(k))'];
Mx3 = 'The solution is:';
Mx4 = 'The error estimate for p is ± ';
clc,echo off, diary output,...
disp(''), disp(Mx1),disp(''), disp(Mx2), disp(points'),...
disp('Iteration converged quadratically to the root.'),...
disp(''),disp(Mx3),disp(''),disp('p = '),...
disp(p2),disp('f(p) = '),disp(y2),...
disp([Mx4,num2str(err)]),diary off,echo on

pause % Press any key to perform Muller's iteration.

clc;
% Place the starting values in p0 p1 and p2

% Place the abscissa tolerance in delta

% Place the ordinate tolerance in epsilon

% Place the number of iterations in max1

p0 = 1.4;
p1 = 1.3;
p2 = 1.2;
delta = 1e-12;
epsilon = 1e-12;
max1 = 10;

[p2,y2,err,P] = muller('f',p0,p1,p2,delta,epsilon,max1);

pause % Press any key for the list of iterations.

clc; clg;
a = 0.975;
b = 1.425;
c = -0.05;
d = 0.55;
h = (b-a)/100;
X = a:h:b;
Y = f(X);
max1 = length(P);
n0 = min(7,max1);
X0 = P(1:n0);
Z0 = zeros(1,n0);
axis([a b c d]);...
plot(X,Y,'-g',X0,Z0,'or');...
hold on;...
plot([a b],[0 0],'b',[0 0],[c d],'b');...
xlabel('x');...
ylabel('y');...
title('Graphical analysis for Muller`s method.');...
grid;...
axis;...
hold off;...
shg; pause % Press any key to continue.

J = 1:max1;
Yp = f(P);
points = [J;P;Yp];
clc,echo off, diary on,...
disp(''), disp(Mx1),disp(''), disp(Mx2), disp(points'),...
disp('Iteration converged quadratically to the root.'),...
disp(''),disp(Mx3),disp(''),disp('p = '),...
disp(p2),disp('f(p) = '),disp(y2),...
disp([Mx4,num2str(err)]),diary off,echo on