diff --git a/examples/maxwell/weak_FEM_BEM_coupling/README.md b/examples/maxwell/weak_FEM_BEM_coupling/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..5e876f9403adc096aeecae03e54f35633d0af697
--- /dev/null
+++ b/examples/maxwell/weak_FEM_BEM_coupling/README.md
@@ -0,0 +1,61 @@
+
+#The electromagnetic transmission-scattering problem
+
+Dependencies:
+* GmshFEM, commit: bd0fb7ebb6b76f0541e25212b9cbcc42e7a0d0b7
+* Gmsh, commit: 2db112186d43fab77fca9d2b546c8f76b3ea6169
+* GmshDdm, commit: 1e539b2a25293011fd80a8b0600a2a55c27fdfbf
+* Bempp-cl, commit: 293135531f385f39ba29cf1c52cc2dba84614e6b
+
+## Method overview:
+
+The method is an efficient weak coupling formulation between the boundary element method
+and the high-order finite element method. The approach is based on the use of a non-overlapping
+domain decomposition method involving quasi-optimal transmission operators. The associated
+transmission boundary conditions are constructed through a localization process based on complex
+rational Padé approximants of the nonlocal Magnetic-to-Electric operators (see I. Badia, B. Caudron, X. Antoine, C. Geuzaine. "A well-conditioned weak coupling of boundary element and high-order finite element methods for time-harmonic electromagnetic scattering by inhomogeneous objects". SIAM Journal on Scientific Computing, Society for Industrial and Applied Mathematics, 2022)
+
+## Problem description
+Let ```math \Omega_{-} ``` be a unit sphere (the scatter) and ```math \Omega_{+} ``` the exterior domain.
+We consider an incident electromagnetic plane wave propagating in ```math \Omega_{+} ```. The total wave field ```math \mathbf{E}``` verifies the following  three-dimensional electromagnetic transmission-scattering problem:
+```math
+\begin{subequations}
+\begin{align}
+\curl \left(\frac{1}{k_{-}\mathcal{Z}_{-}}\curl\E\right) -  k_{-}\mathcal{Z}_{-}^{-1}\E &=  \textbf{0}, \text{ ~in~ } \Omega_{-},  \\
+\curl\curl\E - k_{+}^{2}\E & =  \textbf{0},  \text{ ~in~ } \Omega_{+}, \\
+\frac{1}{\iota k_{+}}\curl \textbf{E}_{\textbf{sc}} \times \frac{\textbf{x}}{\left\lVert\textbf{x} \right\rVert}-\textbf{E}_{\textbf{sc}} &=   \underset{r \to +\infty}{\mathcal{O}}(r^{-2}), \\
+\boldsymbol{\gamma}^{-}_{t}\E &=\boldsymbol{\gamma}^{+}_ {t}\E, \text{ ~on~ } \Gamma, \\
+\boldsymbol{\gamma}^{-}_{t} \left(\frac{1}{k_{-}\mathcal{Z}_{-}}\curl\E\right) &=\frac{1}{k_{+}\mathcal{Z}_{+}}\boldsymbol{\gamma}^{+}_ {t}\left(\curl\E\right), \text{ ~on~ } \Gamma,
+\end{align}
+\end{subequations}
+```
+with ```math k_\pm``` and ```math \mathcal{Z}_\pm``` the wavenumbers and the impedances associated with  ```math \Omega_\pm``` respectively.
+
+## Installation
+
+Run
+
+```
+  mkdir build && cd build
+  cmake ..
+  make
+```
+
+
+## Usage
+
+Simply run:
+```
+  ./example [PARAM]
+```
+with `[PARAM]`:
+* `-pade [x]`  (with `[x]` equal to 1 or 0)  is a parameter to indicate to use Padé or not
+* `-pointsByWl [x]` where `[x]` is the mesh size.
+* `-FEorderAlpha [x]` is the order of the finite element hierarchical basis of the interior electric field.
+* `-FEorderAlpha [x]` is the order of the finite element hierarchical basis of some fields.
+* `-help` to display all other available parameters
+
+The folder contains
+* the main script (`main.cpp`)
+* A python script for the BEM part (`bem.py`)
+* The pybind11 folder for information exchange between python and C++