function p6
  disp('----------------------------------------')
  disp('Known multiplicity:')
  disp('----------------------------------------')
  newton_times_m(1, @f, @fprime, 3, 0, 3);
  disp('----------------------------------------')
  disp('Known second derivative:')
  disp('----------------------------------------')
  newton_second_deriv(1, @f, @fprime, @f2prime, 3, 0);
end

function y = f(x)
  y = x*(1-cos(x));
end
function y = fprime(x)
  y = x*sin(x)-cos(x)+1;
end
function y = f2prime(x)
  y = 2*sin(x)+x*cos(x);
end

function newton_times_m(x, f, fprime, itcount, true, mult)
  disp(sprintf('Iteration:%d value:%.10g error:%.10e', 0, x, abs(x-true)))
  for i=1:itcount
    x = x - mult*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 newton_second_deriv(x, f, fprime, f2prime, 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)/(fprime(x)^2-f(x)*f2prime(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
