diff --git a/ufl/algorithms/__init__.py b/ufl/algorithms/__init__.py
index 2e7fc3a50..86d55c4bc 100644
--- a/ufl/algorithms/__init__.py
+++ b/ufl/algorithms/__init__.py
@@ -36,7 +36,6 @@
     "expand_derivatives",
     "extract_coefficients",
     "extract_base_form_operators",
-    "extract_external_operators",
     "strip_variables",
     "strip_terminal_data",
     "replace_terminal_data",
@@ -73,7 +72,6 @@
     extract_arguments,
     extract_coefficients,
     # extract_arguments_and_coefficients,
-    extract_external_operators,
     extract_base_form_operators,
     extract_elements,
     extract_unique_elements,
diff --git a/ufl/algorithms/analysis.py b/ufl/algorithms/analysis.py
index fc34984de..9e704e8b3 100644
--- a/ufl/algorithms/analysis.py
+++ b/ufl/algorithms/analysis.py
@@ -144,12 +144,6 @@ def extract_coefficients(a):
     return sorted_by_count(extract_type(a, BaseCoefficient))
 
 
-def extract_external_operators(a):
-    """Build a sorted list of all external operators in a,
-    which can be a Form, Integral or Expr."""
-    return sorted_by_count(extract_type(a, ExternalOperator))
-
-
 def extract_constants(a):
     """Build a sorted list of all constants in a"""
     return sorted_by_count(extract_type(a, Constant))
diff --git a/ufl/algorithms/apply_function_pullbacks.py b/ufl/algorithms/apply_function_pullbacks.py
index f7091d497..b7790c7fb 100644
--- a/ufl/algorithms/apply_function_pullbacks.py
+++ b/ufl/algorithms/apply_function_pullbacks.py
@@ -158,11 +158,6 @@ def __init__(self):
     def terminal(self, t):
         return t
 
-    def external_operator(self, e):
-        c = e.result_coefficient(unpack_reference=False)
-        result_coefficient = apply_single_function_pullbacks(c, e.ufl_element())
-        return e._ufl_expr_reconstruct_(*e.ufl_operands, result_coefficient=result_coefficient)
-
     @memoized_handler
     def form_argument(self, o):
         # Represent 0-derivatives of form arguments on reference
diff --git a/ufl/algorithms/check_arities.py b/ufl/algorithms/check_arities.py
index db62735a0..3f2f64a7d 100644
--- a/ufl/algorithms/check_arities.py
+++ b/ufl/algorithms/check_arities.py
@@ -42,10 +42,6 @@ def nonlinear_operator(self, o):
 
     expr = nonlinear_operator
 
-    def external_operator(self, o):
-        # FIXME: Uses map_expr_dag instead of self.arguments because of the Conj case
-        return tuple(map_expr_dag(self, arg, compress=False)[0] for arg in o.argument_slots(outer_form=True))
-
     def sum(self, o, a, b):
         if a != b:
             raise ArityMismatch(f"Adding expressions with non-matching form arguments {_afmt(a)} vs {_afmt(b)}.")
diff --git a/ufl/algorithms/estimate_degrees.py b/ufl/algorithms/estimate_degrees.py
index e63eeb9e1..53dae9354 100644
--- a/ufl/algorithms/estimate_degrees.py
+++ b/ufl/algorithms/estimate_degrees.py
@@ -74,9 +74,6 @@ def coefficient(self, v):
             d = self.default_degree
         return d
 
-    def external_operator(self, v):
-        return self.coefficient(v.result_coefficient())
-
     def _reduce_degree(self, v, f):
         """Reduces the estimated degree by one; used when derivatives
         are taken. Does not reduce the degree when TensorProduct elements
diff --git a/ufl/algorithms/replace.py b/ufl/algorithms/replace.py
index f6763c89b..85970e5d0 100644
--- a/ufl/algorithms/replace.py
+++ b/ufl/algorithms/replace.py
@@ -9,8 +9,7 @@
 #
 # Modified by Anders Logg, 2009-2010
 
-from ufl.core.external_operator import ExternalOperator
-from ufl.classes import CoefficientDerivative, Interp, Form
+from ufl.classes import CoefficientDerivative, Interp, ExternalOperator, Form
 from ufl.constantvalue import as_ufl
 from ufl.corealg.multifunction import MultiFunction
 from ufl.algorithms.map_integrands import map_integrand_dags
@@ -37,29 +36,6 @@ def ufl_type(self, o, *args):
         except KeyError:
             return self.reuse_if_untouched(o, *args)
 
-    """
-    def external_operator(self, o):
-        try:
-            o = self.mapping[o]
-            coeff = o.result_coefficient(unpack_reference=False)
-        except KeyError:
-            coeff = replace(o.result_coefficient(unpack_reference=False), self.mapping)
-        except AttributeError:
-            # ExternalOperator is replaced by something that is not an ExternalOperator
-            return o
-        new_ops = tuple(replace(op, self.mapping) for op in o.ufl_operands)
-        # Fix this
-        if type(new_ops[0]).__name__ == 'Coefficient' and type(o.ufl_operands[0]).__name__ == 'Function':
-            new_ops = o.ufl_operands
-
-        # Does not need to use replace on the 0-slot argument of external operators (v*) since
-        # this can only be a Coargument or a Cofunction, so we directly check into the mapping.
-        # Also, replace is built to apply on Expr, Coargument is not an Expr.
-        new_args = tuple(replace(arg, self.mapping) if not isinstance(arg, Coargument) else self.mapping.get(arg, arg)
-                         for arg in o.argument_slots())
-        return o._ufl_expr_reconstruct_(*new_ops, result_coefficient=coeff, argument_slots=new_args)
-    """
-
     def external_operator(self, o):
         o = self.mapping.get(o) or o
         if isinstance(o, ExternalOperator):
diff --git a/ufl/form.py b/ufl/form.py
index 4707d7b4e..f8b84c227 100644
--- a/ufl/form.py
+++ b/ufl/form.py
@@ -10,7 +10,6 @@
 # Modified by Anders Logg, 2009-2011.
 # Modified by Massimiliano Leoni, 2016.
 # Modified by Cecile Daversin-Catty, 2018.
-# Modified by Nacime Bouziani, 2020.
 
 import warnings
 from collections import defaultdict
@@ -194,6 +193,7 @@ def __call__(self, *args, **kwargs):
         Replaces form.arguments() with given positional arguments in
         same number and ordering. Number of positional arguments must
         be 0 or equal to the number of Arguments in the form.
+
         The optional keyword argument coefficients can be set to a dict
         to replace Coefficients with expressions of matching shapes.