diff --git a/coremltools/converters/mil/frontend/tensorflow/ops.py b/coremltools/converters/mil/frontend/tensorflow/ops.py
index 1f45e821e..ebb07e0ff 100644
--- a/coremltools/converters/mil/frontend/tensorflow/ops.py
+++ b/coremltools/converters/mil/frontend/tensorflow/ops.py
@@ -2146,7 +2146,14 @@ def Sigmoid(context, node):
 @register_tf_op
 def Softplus(context, node):
     x = context[node.inputs[0]]
-    x = mb.softplus(x=x, name=node.name)
+    # Numerically stable softplus: max(x, 0) + log(1 + exp(-|x|)).
+    # Since -|x| <= 0, exp(-|x|) is always in (0, 1], so no overflow occurs
+    # in any precision.  This avoids fp16 overflow on the Apple Neural Engine
+    # at x > ~10.4 (see issues #2687, #2747).
+    abs_x = mb.abs(x=x)
+    exp_val = mb.exp(x=mb.mul(x=-1.0, y=abs_x))
+    log_val = mb.log(x=mb.add(x=1.0, y=exp_val))
+    x = mb.add(x=mb.maximum(x=x, y=0.0), y=log_val, name=node.name)
     context.add(node.name, x)
 
 
diff --git a/coremltools/converters/mil/frontend/tensorflow/test/test_ops.py b/coremltools/converters/mil/frontend/tensorflow/test/test_ops.py
index af78a8475..9fcfd52be 100644
--- a/coremltools/converters/mil/frontend/tensorflow/test/test_ops.py
+++ b/coremltools/converters/mil/frontend/tensorflow/test/test_ops.py
@@ -432,6 +432,92 @@ def build_model(x):
             backend=backend,
         )
 
+    @pytest.mark.parametrize(
+        "compute_unit, backend",
+        itertools.product(
+            [ct.ComputeUnit.CPU_ONLY, ct.ComputeUnit.ALL],
+            backends,
+        ),
+    )
+    def test_softplus_large_values(self, compute_unit, backend):
+        """Verify softplus produces correct output for large inputs (#2747)."""
+        input_shape = (1, 6)
+
+        @make_tf_graph([input_shape])
+        def build_model(x):
+            return tf.math.softplus(x)
+
+        model, inputs, outputs = build_model
+        input_values = [np.array([[15.0, 20.0, -5.0, 0.0, 1.0, 10.5]], dtype=np.float32)]
+        input_dict = dict(zip(inputs, input_values))
+
+        self.run_compare_tf(
+            model,
+            input_dict,
+            outputs,
+            compute_unit=compute_unit,
+            backend=backend,
+        )
+
+    def test_softplus_stable_decomposition(self):
+        """Verify softplus is decomposed into overflow-safe ops for fp16 (#2747).
+
+        The native ``mb.softplus`` op overflows in fp16 on the Apple Neural Engine
+        for x > ~10.4 (exp(x) exceeds 65504).  The stable decomposition
+        ``max(x, 0) + log(1 + exp(-|x|))`` avoids this because exp(-|x|) is
+        always in (0, 1].
+
+        On CPU the native op produces correct output regardless of this fix
+        (the CPU runtime computes in fp32 internally even when the model uses
+        fp16 compute precision — verified empirically).  The overflow only
+        manifests on the ANE, which CI cannot exercise.  We therefore:
+
+        1. Demonstrate the fp16 overflow mathematically with numpy (the same
+           arithmetic the ANE uses), proving the naive formula is unsafe.
+        2. Assert the converter emits the stable decomposition (no native
+           ``softplus`` ops, at least one ``exp`` op), which is the same
+           approach used in PR #2725 for the PyTorch frontend.
+        """
+        # --- 1. Numpy proof: naive fp16 softplus overflows for x > ~10.4 ---
+        x_val = np.float16(15.0)
+        naive = np.log(np.float16(1.0) + np.exp(x_val))
+        assert not np.isfinite(naive), (
+            f"Naive fp16 softplus(15) should overflow to inf, got {naive}"
+        )
+        stable = np.maximum(x_val, np.float16(0.0)) + np.log(
+            np.float16(1.0) + np.exp(-np.abs(x_val))
+        )
+        assert np.isfinite(stable) and stable == np.float16(15.0), (
+            f"Stable fp16 softplus(15) should be 15.0, got {stable}"
+        )
+
+        # --- 2. Conversion-only graph assertion ---
+        input_shape = (1, 6)
+
+        @make_tf_graph([input_shape])
+        def build_model(x):
+            return tf.math.softplus(x)
+
+        model, inputs, outputs = build_model
+        mlmodel = ct.convert(
+            model,
+            inputs=[ct.TensorType(name=inputs[0], shape=input_shape, dtype=np.float32)],
+            compute_precision=ct.precision.FLOAT16,
+            convert_to="mlprogram",
+        )
+        prog = mlmodel._mil_program
+        softplus_ops = prog.find_ops(op_type="softplus")
+        assert len(softplus_ops) == 0, (
+            f"Expected no native 'softplus' ops after stable decomposition, "
+            f"but found {len(softplus_ops)}. The converter should replace "
+            f"softplus with max(x,0)+log(1+exp(-|x|)) for fp16 safety."
+        )
+        exp_ops = prog.find_ops(op_type="exp")
+        assert len(exp_ops) >= 1, (
+            "Expected at least one 'exp' op from the stable decomposition, "
+            "but found none."
+        )
+
     @pytest.mark.parametrize(
         "compute_unit, backend, rank_and_axes",
         itertools.product(