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.10 (Newton-Raphson Method in 2-Dimensions).
% Section 2.7, Newton's Method for Systems, Page 116
echo on; clc; format long; hold off; clear
% This program implements the Newton-Raphson method

% for solving 0 = f1(x,y) and 0 = f2(x,y).

% Use the vector notation X = (x y).

% Define and store the functions f1(X) and f2(X) and

% define and store the functions F(X) and J(X) in the

% M-files f1.m f2.m F.m and J.m respectively.

pause % Press any key to continue.

clc;
% function z = f1(x,y)
% z = x.^2 - 2.*x - y + 0.5;

% function z = f2(x,y)
% z = x.^2 + 4.*y.^2 - 4;

% function Z = F(X)
% x = X(1); y = X(2);
% Z = [f1(x,y); f2(x,y)];

% function W = J(X)
% x = X(1); y = X(2);
% W = [(2.*x - 2) (-1);
% (2.*x) (8.*y)];

delete f1.m
diary f1.m; disp('function z = f1(x,y)');...
disp('z = x.^2 - 2.*x - y + 0.5;');...
diary off;

delete f2.m
diary f2.m; disp('function z = f2(x,y)');...
disp('z = x.^2 + 4.*y.^2 - 4;');...
diary off;

delete F.m
diary F.m; disp('function Z = F(X)');...
disp('x = X(1); y = X(2);');...
disp('Z = [f1(x,y); f2(x,y)];');...
diary off;

delete J.m
diary J.m; disp('function W = J(X)');...
disp('x = X(1); y = X(2);');...
disp('W = [(2.*x - 2) (-1);');...
disp(' (2.*x) (8.*y)];');...
diary off;

% Remark. f1.m, f2.m, F.m, J.m and new2dim.m are used for Algorithm 2.10
f1(0,0); f2(0,0); F([0 0]); J([0 0]); % Test for files f1.m f2.m F.m J.m
pause % Press any key to graph 0 = f1(x,y) and 0 = f2(x,y).

clc; clg;
a = -2.0;
b = 2.5;
c = -1.0;
d = 1.5;
hx = (b-a)/31;
hy = (d-c)/31;
[X Y] = meshdom(a:hx:b, c:hy:d);...
Z = f1(X,Y);...
W = f2(X,Y);...
contour(Z,[0 0],a:hx:b,c:hy:d,'-g');...
hold on;...
contour(W,[0 0],a:hx:b,c:hy:d,'-g');...
title('Implicit plot of 0 = f1(x,y) and 0 = f2(x,y).');...
grid;...
hold off;...
shg; pause % Press any key to perform Newton-Raphson iteration.

clc;
% Place the starting value in [p0 q0]

% Place the tolerance for P in delta

% Place the tolerance for F(P) in epsilon

% Place the maximum number of iterations in max1

p0 = 2.0;
q0 = 0.25;
P0 = [p0 q0];
delta = 1e-12;
epsilon = 1e-12;
max1 = 40;

[P0,F0,err,P] = new2dim('F','J',P0,delta,epsilon,max1);

pause % Press any key for the list of iterations.

clg; clc;
a = 1.86;
b = 2.02;
c = 0.24;
d = 0.34;
hx = (b-a)/20;
hy = (d-c)/20;
[X Y] = meshdom(a:hx:b, c:hy:d);...
Z = f1(X,Y);...
W = f2(X,Y);...
contour(Z, [0 0],a:hx:b,c:hy:d,'-g');...
grid;...
hold on;...
contour(W, [0 0],a:hx:b,c:hy:d,'-g');...
X0 = P(:,1);...
Y0 = P(:,2);...
plot(X0,Y0,'or');...
title('Graphical presentation of the Newton-Raphson iteration.');...
hold off;...
shg; pause % Press any key to continue.

Mx1 = 'Iterations for Newton`s method.';
Mx2 = ' p(k) q(k)';
Mx3 = 'The solution is:';
Mx4 = 'The error estimate for P is ';
Mx5 = num2str(err);
Mx6 = [Mx4,'[± ',Mx5,', ± ',Mx5,']'];
clc,echo off, diary output,...
disp(''), disp(Mx1),disp(''), disp(Mx2),disp(P),...
disp('Iteration converged to the root.'),...
disp(''),disp(Mx3),disp(''),disp('P = '),...
disp(P0),disp('F(P) = '),disp(F0),...
disp(''),disp([Mx6]),diary off,echo on

pause % Press any key to perform Newton-Raphson iteration.

clc;
% Place the starting value in [p0 q0]

% Place the tolerance for P in delta

% Place the tolerance for F(P) in epsilon

% Place the maximum number of iterations in max1

p0 = -0.1;
q0 = 0.9;
P0 = [p0 q0];
delta = 1e-12;
epsilon = 1e-12;
max1 = 40;

[P0,F0,err,P] = new2dim('F','J',P0,delta,epsilon,max1);

pause % Press any key for the list of iterations.

clg; clc;
a = -0.26;
b = -0.08;
c = 0.88;
d = 1.02;
hx = (b-a)/20;
hy = (d-c)/20;
[X Y] = meshdom(a:hx:b, c:hy:d);...
Z = f1(X,Y);...
W = f2(X,Y);...
contour(Z,[0 0],a:hx:b,c:hy:d,'-g');...
grid;...
hold on;...
contour(W,[0 0],a:hx:b,c:hy:d,'-g');...
X0 = P(:,1);...
Y0 = P(:,2);...
plot(X0,Y0,'or');...
title('Graphical presentation of the Newton-Raphson iteration.');...
hold off;...
shg; pause % Press any key to continue.

clc,echo off, diary output,...
disp(''), disp(Mx1),disp(''), disp(Mx2),disp(P),...
disp('Iteration converged to the root.'),...
disp(''),disp(Mx3),disp(''),disp('P = '),...
disp(P0),disp('F(P) = '),disp(F0),...
disp(''),disp([Mx6]),diary off,echo on