function p3
  newton(4, @f, @fprime, 3, pi);
  disp(sprintf('For reference: pi=%.10g', pi))
end

function y = f(x)
  y = sin(x);
end
function y = fprime(x)
  y = cos(x);
end

function newton(x, f, fprime, itcount, true)
  disp(sprintf('Iteration:%d value:%.10g error:%.10e', 0, x, abs(x-true)))
  for i=1:itcount
    x = x - f(x)/fprime(x);
    disp(sprintf('Iteration:%d value:%.10g error:%.10e', i, x, abs(x-true)))
    errors(i) = abs(x-true);
  end

  order_of_convergence(errors)
end

function order_of_convergence(seq)
  for at=1:length(seq)-1
    order = log(seq(at+1)) / log(seq(at+0));
    disp(sprintf('Apparent order of convergence at iteration %d:%g', at, order))
  end
end
