f = @(x) x^6 - x - 1;
f_prime = @(x) 6*x^5 - 1;
tol = 1/1000000;
x0 = 2;
expected_root = 1.1347;
unit_test_equals("Newton's method agrees with Haskell", ...
		 expected_root, ...
		 newtons_method(f, f_prime, tol, x0));


f1 = @(u) u(1)^2 + u(1)*u(2)^3 - 9;
f2 = @(u) 3*u(1)^2*u(2) - u(2)^3 - 4;
f = @(u) [f1(u); f2(u)];
## The partials for the Jacobian.
f1x = @(u) 2*u(1) + u(2)^3;
f1y = @(u) 3*u(1)*u(2)^2;
f2x = @(u) 6*u(1)*u(2);
f2y = @(u) 3*u(1)^2 - 3*u(2)^2;
## f_prime == Jacobian.
f_prime = @(u) [ f1x(u), f1y(u); f2x(u), f2y(u) ];
tol = 1 / 10^12;
u0 = [1.2; 2.5];
expected_root = [1.33635; 1.75424];
[actual_root, iterations] = newtons_method(f, f_prime, tol, u0);
unit_test_equals("Homework #3 problem #4 root is correct", ...
		 expected_root, ...
		 actual_root);
