Source code in the pl/i programming language

f: procedure (x) returns (float (18));
   declare x float (18);
   return (x**3 - 3*x**2 + 2*x );
end f;

declare eps float, (x, y) float (18);
declare dx fixed decimal (15,13);

eps = 1e-12;

do dx = -5.03 to 5 by 0.1;
   x = dx;
   if sign(f(x)) ^= sign(f(dx+0.1)) then
      call locate_root;
end;

locate_root: procedure;
   declare (left, mid, right) float (18);

   put skip list ('Looking for root in [' || x, x+0.1 || ']' );
   left = x; right = dx+0.1;
   PUT SKIP LIST (F(LEFT), F(RIGHT) );
   if abs(f(left) ) < eps then
      do; put skip list ('Found a root at x=', left); return; end;
   else if abs(f(right) ) < eps then
      do; put skip list ('Found a root at x=', right); return; end;
   do forever;
      mid = (left+right)/2;
      if sign(f(mid)) = 0 then
         do; put skip list ('Root found at x=', mid); return; end;
      else if sign(f(left)) ^= sign(f(mid)) then
         right = mid;
      else
         left = mid;
      /* put skip list (left || right); */
      if abs(right-left) < eps then
         do; put skip list ('There is a root near ' ||
            (left+right)/2); return;
         end;
   end;
end locate_root;