% Newton-Raphson root finding - single function
%


function func = funct(x)
  func = x^4 - 5*x^3 - 124*x^2 + 380*x + 1200;
end


function deriv = derivfunct(x)
  deriv = 4*x^3 - 15*x^2 - 248*x + 380;
end


% now proceed with finding a root:

x = 1;% Enter initial guess of root location;

itermax = 100;      % max # of iterations
iter = 0;
errmax = 0.001;     % convergence tolerance
error = 1;

write("iter  ","xnew  ","error ");

while error > errmax & iter < itermax

  iter = iter + 1;

  f=funct(x);
  fprime= derivfunct(x);

  if fprime == 0
    write("ERROR: deriv(x) = 0; cant divide by zero")
    break;
  end;

  xnew = x - f / fprime;    % here is new root estimate

  error = abs((xnew - x)/xnew) * 100;  % find change from previous
  write(iter,"  ",xnew,"  ",error);
  x = xnew;      % set up for next iteration

end