e413.md

Problem 4.13

e413.m

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

e413f.m

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