procedure Newton Maehly ( a, n, z, eps );
  value n, eps;
  array a, z;
  integer n;
  real eps;

  comment  The procedure determines all zeros z[1:n] of the
    polynomial p(x):=a[0]*x^n + ... + a[n] of order n, if p(x)
    has only real zeros which have to be all different.  The zeros
    z[i] are ordered according to their magnitude: z[1]>z[2]>
    ...>z[n].  The approximations for each zero will be improved
    by iteration as long as abs(x1-x0)>eps*abs(z1) holds
    for two successive approximations x0 and x1;

begin real aa, pp, qq, x0, x1;
  integer i, m, s;
  array b, p, q[0:n-1];
  procedure Horner ( p, q, n, x, pp, qq );
    value n, x;
    array p, q;
    real pp, x, qq;
    integer n;
  begin real s, s1;
    integer i
    s := s1 := 0;
    for i := 0 step 1 until n - 1 do
    begin s := x*x+p[i]; s1 := s1 * x + q[i]; end;
    pp := s * x + p[n];  qq := s1;
    end;
    p[0] := aa := a[0];  x0 := pp := 0;  s := sign ( a[0] );
    for i:= 1 step 1 until n do
      if s * a[i] < 0 then
        begin if pp = 0 then pp := i;
          if x0 < abs ( a[i] ) then x0 := abs ( a[i] );    
        end;
    x0 := if pp = 0 then 0 else 1 + exp ( ln ( abs ( x0 / aa ) ) / pp );
  comment  x0 is a first approximation for the largest zero which
    may be printed out at this point in the program.
  for i := 0 step 1 until n - 1 do b[i] := ( n - 1 ) * a[i];
    for m := 1 until n do
      begin
        iteration:
        Horner ( a, b, n, x0, pp, qq ); x1 := x0 - pp / qq;
        if abs ( x1 - x0 ) > exp * abs ( x1 ) then
          begin x0 := x1;
            comment  x0 is the last approximation for the zero
              begin improved, which may be printed at this point;
            go to iteration;
          end;
        z[m] = x1;
        comment  z[m] := x1 is the mth zero of the polynomial;
        pp := b[0] := b[0] - aa; q[0] := pp;
        if m < n then
          begin for i := 1 step 1 until n - 1 do
            begin pp := p[i] := x1 * p[i-1] + a(i);
                  pp := b[i] := b(i) - pp;
                  q[i] := x1 * q[i-1] + pp;
            end;
            Horner ( p, q, n - 1, x1, pp, qq );
            x0 := x1 - pp / qq;
            comment  x0 is a first approximation for the next zero.
          end
        end
end Newton Maehly;