abserr = 1e-8;
relerr = 1e-6;
figure;
fplot(@(x) cot(x) - (x^2 - 1) / (2*x), [0 10, -10, 10]);
for n = 1:3
b = (n - 1) * pi + 1/(2*n);
c = (n - 1) * pi + pi / 2;
disp([b c]);
[b,c,residual,flag] = Zero('e413f',b,c,abserr,relerr);
% Check flag and print results.
fprintf('flag = %i \n',flag)
if flag == 0
fprintf('Computed a root b = %e \n',b);
elseif flag == 1
fprintf('Too much work\n');
fprintf('There is a root in [b,c] with \n');
fprintf('b = %e, c = %e \n',b,c);
elseif flag == 2
fprintf('Computed a pole b = %e \n',b);
end
fprintf('The residual f(b) = %e \n',residual);
end
function y = e413f( x )
y = cot(x) - (x^2 - 1) / (2*x);
return
stdout:
e413
flag = 0
Computed a root b = 1.306543e+00
The residual f(b) = -3.140715e-07
3.3916 4.7124
flag = 0
Computed a root b = 3.673194e+00
The residual f(b) = 1.567060e-06
6.4499 7.8540
flag = 0
Computed a root b = 6.584620e+00
The residual f(b) = 1.362247e-06
9.5498 10.9956
flag = 0
Computed a root b = 9.631683e+00
The residual f(b) = 3.198718e-05