diff --git a/.gitignore b/.gitignore
index ff2274027fbec47efc58c49bc2ce58f66d75552a..82c341ae600d2511ffde37d4f316e1dab5d353dc 100644
--- a/.gitignore
+++ b/.gitignore
@@ -20,6 +20,7 @@ examples/*.pdf
 
 *.vtu
 *.vts
+*.dat
 
 .cache
 
diff --git a/experiments/cahn-hilliard.py b/experiments/cahn-hilliard.py
deleted file mode 100644
index 8393f744caae96c798dd293be299b16127df1dc4..0000000000000000000000000000000000000000
--- a/experiments/cahn-hilliard.py
+++ /dev/null
@@ -1,132 +0,0 @@
-import numpy as np
-import pyopencl as cl
-import pyopencl.clmath  # noqa
-
-from meshmode.discretization import Discretization
-from meshmode.discretization.poly_element import \
-        InterpolatoryQuadratureSimplexGroupFactory
-
-from pytential import bind, sym, norm  # noqa
-from pytential.target import PointsTarget
-
-# {{{ set some constants for use below
-
-nelements = 20
-bdry_quad_order = 4
-mesh_order = bdry_quad_order
-qbx_order = bdry_quad_order
-bdry_ovsmp_quad_order = 4*bdry_quad_order
-fmm_order = 8
-
-# }}}
-
-
-def main():
-    import logging
-    logging.basicConfig(level=logging.INFO)
-
-    cl_ctx = cl.create_some_context()
-    queue = cl.CommandQueue(cl_ctx)
-
-    from meshmode.mesh.generation import ellipse, make_curve_mesh
-    from functools import partial
-
-    mesh = make_curve_mesh(
-                partial(ellipse, 2),
-                np.linspace(0, 1, nelements+1),
-                mesh_order)
-
-    pre_density_discr = Discretization(
-            cl_ctx, mesh,
-            InterpolatoryQuadratureSimplexGroupFactory(bdry_quad_order))
-
-    from pytential.qbx import (
-            QBXLayerPotentialSource, QBXTargetAssociationFailedException)
-    qbx, _ = QBXLayerPotentialSource(
-            pre_density_discr, fine_order=bdry_ovsmp_quad_order, qbx_order=qbx_order,
-            fmm_order=fmm_order,
-            expansion_disks_in_tree_have_extent=True,
-            ).with_refinement()
-    density_discr = qbx.density_discr
-
-    from pytential.symbolic.pde.cahn_hilliard import CahnHilliardOperator
-    chop = CahnHilliardOperator(
-            # FIXME: Constants?
-            lambda1=1.5,
-            lambda2=1.25,
-            c=1)
-
-    unk = chop.make_unknown("sigma")
-    bound_op = bind(qbx, chop.operator(unk))
-
-    # {{{ fix rhs and solve
-
-    nodes = density_discr.nodes().with_queue(queue)
-
-    def g(xvec):
-        x, y = xvec
-        return cl.clmath.atan2(y, x)
-
-    bc = sym.make_obj_array([
-        # FIXME: Realistic BC
-        g(nodes),
-        -g(nodes),
-        ])
-
-    from pytential.solve import gmres
-    gmres_result = gmres(
-            bound_op.scipy_op(queue, "sigma", dtype=np.complex128),
-            bc, tol=1e-8, progress=True,
-            stall_iterations=0,
-            hard_failure=True)
-
-    # }}}
-
-    # {{{ postprocess/visualize
-
-    sigma = gmres_result.solution
-
-    from sumpy.visualization import FieldPlotter
-    fplot = FieldPlotter(np.zeros(2), extent=5, npoints=500)
-
-    targets = cl.array.to_device(queue, fplot.points)
-
-    qbx_stick_out = qbx.copy(target_association_tolerance=0.05)
-
-    indicator_qbx = qbx_stick_out.copy(qbx_order=2)
-
-    from sumpy.kernel import LaplaceKernel
-    ones_density = density_discr.zeros(queue)
-    ones_density.fill(1)
-    indicator = bind(
-            (indicator_qbx, PointsTarget(targets)),
-            sym.D(LaplaceKernel(2), sym.var("sigma")))(
-            queue, sigma=ones_density).get()
-
-    try:
-        fld_in_vol = bind(
-                (qbx_stick_out, PointsTarget(targets)),
-                chop.representation(unk))(queue, sigma=sigma).get()
-    except QBXTargetAssociationFailedException as e:
-        fplot.write_vtk_file(
-                "failed-targets.vts",
-                [
-                    ("failed", e.failed_target_flags.get(queue))
-                    ]
-                )
-        raise
-
-    #fplot.show_scalar_in_mayavi(fld_in_vol.real, max_val=5)
-    fplot.write_vtk_file(
-            "potential.vts",
-            [
-                ("potential", fld_in_vol),
-                ("indicator", indicator),
-                ]
-            )
-
-    # }}}
-
-
-if __name__ == "__main__":
-    main()
diff --git a/pytential/symbolic/pde/cahn_hilliard.py b/pytential/symbolic/pde/cahn_hilliard.py
index 3263f18e1b9c144a743335107339770a12425c30..00d17038588beeb65ea1e29eeb0cf1ca55207e8c 100644
--- a/pytential/symbolic/pde/cahn_hilliard.py
+++ b/pytential/symbolic/pde/cahn_hilliard.py
@@ -34,10 +34,60 @@ from functools import partial
 
 
 class CahnHilliardOperator(L2WeightedPDEOperator):
-    def __init__(self, lambda1, lambda2, c):
-        self.lambdas = (lambda1, lambda2)
+    def __init__(self, b, c):
+        self.b = b
         self.c = c
 
+        # Issue:
+        # - let crat = np.abs(4.*c) / ( b**2 + 1e-12 )
+        # - when crat is close to zero, sqrt(b**2-4*c) is close to abs(b),
+        #   then for b>=0, sqrt(b**2-4*c) - b is inaccurate.
+        # - similarly, when crat is close to one, sqrt(b**2-4*c) is close to zero,
+        #   then for b>0, sqrt(b**2-4*c) + b is inaccurate.
+        # Solution:
+        # - find a criteria for crat to choose from formulae, or
+        # - use the computed root with smaller residual
+        def quadratic_formula_1(a, b, c):
+            return (-b + np.sqrt(b**2 - 4*a*c)) / (2*a)
+
+        def quadratic_formula_2(a, b, c):
+            return (-b - np.sqrt(b**2 - 4*a*c)) / (2*a)
+
+        def citardauq_formula_1(a, b, c):
+            return 2*c / (-b - np.sqrt(b**2-4*a*c))
+
+        def citardauq_formula_2(a, b, c):
+            return 2*c / (-b + np.sqrt(b**2-4*a*c))
+
+        def f(x):
+            return x**2 - b*x + c
+
+        root11 = quadratic_formula_1(1, -b, c)
+        root12 = citardauq_formula_1(1, -b, c)
+        if np.abs(f(root11)) < np.abs(f(root12)):
+            lam1 = np.sqrt(root11)
+        else:
+            lam1 = np.sqrt(root12)
+
+        root21 = quadratic_formula_2(1, -b, c)
+        root22 = citardauq_formula_2(1, -b, c)
+        if np.abs(f(root21)) < np.abs(f(root22)):
+            lam1 = np.sqrt(root21)
+        else:
+            lam2 = np.sqrt(root22)
+
+        if not np.abs(f(lam1**2)) < 1e-12:
+            import warnings
+            warnings.warn("Root lam1's residual is less than 1e-12.")
+            print(np.abs(f(lam1**2)))
+        if not np.abs(f(lam2**2)) < 1e-12:
+            import warnings
+            warnings.warn("Root lam2's residual is less than 1e-12.")
+            print(np.abs(f(lam2**2)))
+
+        # self.lambdas = sorted([lam1, lam2], key=abs, reverse=True)  # biggest first
+        self.lambdas = sorted([lam1, lam2], key=abs, reverse=False)  # smallest first
+
     def make_unknown(self, name):
         return sym.make_sym_vector(name, 2)
 
@@ -45,60 +95,122 @@ class CahnHilliardOperator(L2WeightedPDEOperator):
         if op_map is None:
             op_map = lambda x: x  # noqa: E731
 
-        from sumpy.kernel import HelmholtzKernel
-        hhk = HelmholtzKernel(2, allow_evanescent=True)
-        hhk_scaling = 1j/4
+        from sumpy.kernel import YukawaKernel
+        knl = YukawaKernel(2)
 
         if i == 0:
             lam1, lam2 = self.lambdas
             return (
-                    # FIXME: Verify scaling
-                    -1/(2*np.pi*(lam1**2-lam2**2)) / hhk_scaling
+                    1/(lam1**2-lam2**2)
                     * (
-                        op_map(sym.S(hhk, density, k=1j*lam1,
-                            qbx_forced_limit=qbx_forced_limit))
-                        - op_map(sym.S(hhk, density, k=1j*lam2,
-                            qbx_forced_limit=qbx_forced_limit))))
+                        op_map(sym.S(knl, density, lam=lam1,
+                               qbx_forced_limit=qbx_forced_limit))
+                        - op_map(sym.S(knl, density, lam=lam2,
+                                 qbx_forced_limit=qbx_forced_limit))))
         else:
             return (
-                    # FIXME: Verify scaling
-
-                    -1/(2*np.pi) / hhk_scaling
-                    * op_map(sym.S(hhk, density, k=1j*self.lambdas[i-1],
+                    op_map(sym.S(knl, density, lam=self.lambdas[i-1],
                         qbx_forced_limit=qbx_forced_limit)))
 
-    def representation(self, unknown):
+    def debug_representation(self, unknown, qbx_forced_limit=None):
+        """Return some layer potentials used for debugging
+        """
         sig1, sig2 = unknown
-        S_G = partial(self.S_G, qbx_forced_limit=None)  # noqa: N806
+        S_G = partial(self.S_G, qbx_forced_limit=qbx_forced_limit)  # noqa: N806
 
-        return S_G(1, sig1) + S_G(0, sig2)
+        return sym.make_obj_array([
+            S_G(0, sig1),
+            S_G(1, sig1),
+            S_G(2, sig1),
+            ])
 
-    def operator(self, unknown):
+    def representation(self, unknown, qbx_forced_limit=None):
+        """Return (u, v) in a :mod:`numpy` object array.
+        """
         sig1, sig2 = unknown
         lam1, lam2 = self.lambdas
-        S_G = partial(self.S_G, qbx_forced_limit=1)  # noqa: N806
+        S_G = partial(self.S_G, qbx_forced_limit=qbx_forced_limit)  # noqa: N806
+        laplacian = partial(sym.laplace, 2)
+        d_dx = partial(sym.d_dx, 2)
+        d_dy = partial(sym.d_dy, 2)
+
+        def u(op_map=None):
+            return S_G(1, sig1, op_map=op_map) + S_G(0, sig2, op_map=op_map)
+
+        def v(op_map=None):
+            return lam2**2 * S_G(1, sig1, op_map=op_map) \
+                    - S_G(1, sig2, op_map=op_map) \
+                    + lam1**2 * S_G(0, sig2, op_map=op_map)
+
+        # v = - lap_u + b*u
+        return sym.make_obj_array([
+            u(),
+            -u(op_map=laplacian) + self.b*u(),
+            u(op_map=d_dx),
+            u(op_map=d_dy),
+            u(op_map=laplacian),
+            # self.b*u() - v(),
+            v(),
+            v(op_map=d_dx),
+            v(op_map=d_dy)
+            ])
+
+    def operator(self, unknown, rscale=None, bdry_splitting_parameter=0,
+            right_precon=None, right_precon_inv=None,
+            left_precon=None, left_precon_inv=None):
+        """
+        rscale: this API is deprecated, use left_precon instead.
+        left_precon, right_precon: left/right preconditionors.
+        It only accepts diagonal preconditioners for the moment.
+
+        FIXME: Need to use inner product under inv(L*L')
+        """
+        lam1, lam2 = self.lambdas
+        S_G = partial(self.S_G, qbx_forced_limit=-1)  # noqa: N806
+
+        if rscale is not None:
+            raise NotImplementedError("rscale deprecated, use left_precon instead")
+
+        if right_precon is None:
+            right_precon = lambda x: x  # noqa: E731
+            right_precon_inv = right_precon  # noqa: F841
 
+        if left_precon is None:
+            left_precon = lambda x: x  # noqa: E731
+            left_precon_inv = left_precon  # noqa: F841
+
+        sig1, sig2 = right_precon(unknown)
+
+        # c = 1 / (epsilon * delta_t)
         c = self.c
+        # bdry_splitting_parameter = alpha_tild / epsilon
+        bs = bdry_splitting_parameter
 
         def Sn_G(i, density):  # noqa
             return self.S_G(i, density,
                         qbx_forced_limit="avg",
                         op_map=partial(sym.normal_derivative, 2))
 
-        d = sym.make_obj_array([
-            0.5*sig1,
-            0.5*lam2**2*sig1 - 0.5*sig2
+        d_pre = sym.make_obj_array([
+            - 0.5 * sig1,
+            - 0.5 * lam2**2 * sig1 + 0.5 * sig2
             ])
-        a = sym.make_obj_array([
+        d = left_precon(d_pre)
+
+        # The 2nd bc has a nullspace (why?)
+
+        # Relaxational BC
+        a_pre = sym.make_obj_array([
             # A11
-            Sn_G(1, sig1) + c*S_G(1, sig1)
+            Sn_G(1, sig1) + (c-bs)*S_G(1, sig1)
             # A12
-            + Sn_G(0, sig2) + c*S_G(0, sig2),
+            + Sn_G(0, sig2) + (c-bs)*S_G(0, sig2),
 
             # A21
             lam2**2*Sn_G(1, sig1)
             # A22
             - Sn_G(1, sig2) + lam1**2*Sn_G(0, sig2)
             ])
+        a = left_precon(a_pre)
 
         return d+a
diff --git a/pytential/symbolic/pde/scalar.py b/pytential/symbolic/pde/scalar.py
index 882b1b3c53cb53830bb94796b887a71494fdae32..759cad3aaf1090b44ffc1ac92aefa435ebb2ace0 100644
--- a/pytential/symbolic/pde/scalar.py
+++ b/pytential/symbolic/pde/scalar.py
@@ -308,4 +308,104 @@ class NeumannOperator(L2WeightedPDEOperator):
 # }}}
 
 
+# {{{ robin
+
+class RobinOperator(L2WeightedPDEOperator):
+    """IE operator and field representation for solving Robin boundary
+    value problems with scalar kernels a * u + b * u_n = g (e.g.
+    :class:`sumpy.kernel.LaplaceKernel`)
+
+    .. automethod:: is_unique_only_up_to_constant
+    .. automethod:: representation
+    .. automethod:: operator
+    """
+
+    def __init__(self, kernel,
+            loc_sign,
+            a=1., b=1.,
+            laplace_kernel=0, use_l2_weighting=False,
+            kernel_arguments=None):
+        """
+        :arg a: coefficient for the function value u.
+        :arg b: coefficient for the normal derivative u_n.
+        :arg loc_sign: +1 for exterior, -1 for interior
+        """
+
+        assert a**2 + b**2 > 0
+        assert loc_sign in [-1, 1]
+
+        assert isinstance(kernel, Kernel)
+
+        if kernel_arguments is None:
+            kernel_arguments = {}
+        self.kernel_arguments = kernel_arguments
+
+        self.a = a
+        self.b = b
+        self.loc_sign = loc_sign
+
+        L2WeightedPDEOperator.__init__(self, kernel, use_l2_weighting)
+
+    def is_unique_only_up_to_constant(self):
+        from sumpy.kernel import LaplaceKernel
+        return (
+                isinstance(self.kernel, LaplaceKernel)
+                and self.loc_sign < 0
+                and self.a == 0
+                )
+
+    def representation(self, u, map_potentials=None, qbx_forced_limit=None,
+            **kwargs):
+        sqrt_w = self.get_sqrt_weight()
+        inv_sqrt_w_u = cse(u/sqrt_w)
+
+        if map_potentials is None:
+            def map_potentials(x):  # pylint:disable=function-redefined
+                return x
+
+        kwargs["qbx_forced_limit"] = qbx_forced_limit
+        kwargs["kernel_arguments"] = self.kernel_arguments
+
+        return (1. / (self.a**2 + self.b**2)) * (
+                self.b * map_potentials(
+                    sym.S(self.kernel, inv_sqrt_w_u, **kwargs))
+                -
+                self.a * map_potentials(
+                    sym.D(self.kernel, inv_sqrt_w_u, **kwargs))
+                )
+
+    def operator(self, u):
+        sqrt_w = self.get_sqrt_weight()
+        inv_sqrt_w_u = cse(u/sqrt_w)
+
+        knl = self.kernel
+
+        knl_kwargs = {}
+        knl_kwargs["kernel_arguments"] = self.kernel_arguments
+
+        if self.is_unique_only_up_to_constant():
+            # The interior Neumann operator in this representation
+            # has a nullspace. The mean of the density must be matched
+            # to the desired solution separately. As is, this operator
+            # returns a mean that is not well-specified.
+
+            amb_dim = self.kernel.dim
+            ones_contribution = (
+                    sym.Ones() * sym.mean(amb_dim, amb_dim-1, inv_sqrt_w_u))
+        else:
+            ones_contribution = 0
+
+        return (-self.loc_sign * 0.5 * u
+                + (1. / (self.a**2 + self.b**2)) * sqrt_w * (
+                    self.a * self.b * sym.S(self.kernel, inv_sqrt_w_u, **kwargs)
+                    - self.a**2 * sym.D(self.kernel, inv_sqrt_w_u, **kwargs)
+                    + self.b**2 * sym.Sp(
+                        knl, inv_sqrt_w_u, qbx_forced_limit="avg", **knl_kwargs)
+                    - self.b * self.a * sym.Dp(
+                        knl, inv_sqrt_w_u, qbx_forced_limit="avg", **knl_kwargs)
+                    + ones_contribution
+                    ))
+
+# }}}
+
 # vim: foldmethod=marker