I'm solving an eigenvalue problem with special b.c. mentioned in the paper [*Edge states in gated bilayer-monolayer graphene ribbons and bilayer domain walls*](https://doi.org/10.1063/1.5025937). (See Eq.(3)(6)(7) and the paragraph after (9). ) The paper author claims the problem is solved with COMSOL MULTIPHYSICS. I tried with `NDEigensystem` but failed.

The system involves four coupled ODE [$A1(x),B1(x),A2(x),B2(x)$] defined in the interval $x=\{-\infty,0\}$ and 2 ODE [$A(x),B(x)$] Eqs for $x=\{0,W\}$, The b.c.s are:

$A1(0)=A(0)$  

$B1(0)=B(0)$  

$B2(0)=0$  

$B(W)=0$

and all functions $(A1,B1,A2,B2)$ in the left interval go to zero as $x$ goes to $-\infty$. Here is my trivial try

```
systm1[k_]:=Module[{a=400,hv = 0.65875 10^3, W = 23.86},{vals,funs}=NDEigensystem[{

hv(-B1'[x]-k B1[x]) +50 A1[x],
hv(A1'[x]-k A1[x]) +50 B1[x]+a A2[x],
hv(B2'[x]+k B2[x]) -50 A2[x]+a B1[x],
hv(-A2'[x]+k A2[x]) -50 B2[x],

hv(-B'[x]-k B[x])+50A[x] ,
hv(A'[x]-k A[x])+50B[x] ,

DirichletCondition[{A1[0]==A[0],B1[0]==B[0],B2[0]==0,B[W]==0},True]},{A1[x],B1[x],A2[x],B2[x],A[x],B[x]},{x,-10W,10W},6,Method->{"SpatialDiscretization"->{"FiniteElement",{"MeshOptions"->{MaxCellMeasure->0.2}}}}];
Return[vals];]   

systm1[0.]

```

> 
> NDEigensystem::femcscd: The PDE is convection dominated and the result
> may not be stable. Adding artificial diffusion may help.
> 
> 
> NDEigensystem::fembdcc: Cross-coupling of dependent variables in
> `DirichletCondition[A1 == A, True]` is not supported in this version.
> 
> 
> 

Kindly how can I implement the desired BCs?

---

Update 1
--------

According to the answer by @xzczd, the problem can be reduced to 4 coupled odes considering the continuity at `x=0` so the system with new BCs can be seen as in this picture

[![enter image description here](https://i.stack.imgur.com/0qWoa.png)](https://i.stack.imgur.com/0qWoa.png)

---

Update 2
--------

So, following the answer by @xzczd we consider only the first 4 ode in the question and restrict the `A2` and `B2` (using `If`)to be only in the interval $(-\infty,0)$ which becomes

```
{a = 400, hv = 0.65875 10^3, W = 23.86, k = 0., inf = 10 W, bR = W};
<< NDSolve`FEM`;
meshline = 
  ToElementMesh[ImplicitRegion[-inf <= x <= bR, x], {{-inf, bR}}, 
   IncludePoints -> {{0}}, MaxCellMeasure -> 1];

{val, vec} = NDEigensystem[Flatten@{

     hv (-B1'[x] - k B1[x]) + 50 A1[x],
     hv (A1'[x] - k A1[x]) + 50 B1[x] + If[x < 0, a, 0] A2[x],
     hv (B2'[x] + k B2[x]) - If[x < 0, 50, 0] A2[x] + a B1[x],
     hv (-A2'[x] + k A2[x]) - If[x < 0, 50, 0] B2[x],

     DirichletCondition[B2[x] == 0, x == 0], 
     DirichletCondition[B1[x] == 0, x == bR], 
     DirichletCondition[#, x == -inf] & /@ {A1[x] == 0, B1[x] == 0, 
       A2[x] == 0, B2[x] == 0}}, {A1[x], B1[x], A2[x], 
    B2[x]}, {x} \[Element] meshline, 6];      

```

but the eigenfunctions do not satisfy the some of BCs.

```
func = {A1, B1, A2, B2};
Plot[Abs[vec[[1, #]]], {x, -3 W, W}, PlotRange -> All, Frame -> True, 
   PlotStyle -> Red, PlotLabel -> ToString[func[[#]]]] & /@ Range[4]    

```

[![enter image description here](https://i.stack.imgur.com/prBmZ.png)](https://i.stack.imgur.com/prBmZ.png)

---

Update 3: Exact solution for k=0
--------------------------------

For `k=0`, we can find an analytic solution, I will use here dimensionless units such that `hv=1;a=1`. So for $x<0$ we have these coupled equations

```
-B1'[x] +u A1[x]==En A1[x];
 A1'[x] +u B1[x]+A2[x]==En B1[x];
 B2'[x] -u A2[x]+B1[x]==En A2[x];
-A2'[x] -u B2[x]==En B2[x]; 

```

`En` is the eigenvalue, these equations can be decoupled and we obtain:

```
A1[x_]=(b1 E^(g1 x) g1)/(-En + u) + (b3 E^(g2 x) g2)/(-En + u);
B1[x_]=b1 E^(g1 x) + b3 E^(g2 x);
A2[x_]=(b1 E^(g1 x) (g1^2 + (En - u)^2) + 
 b3 E^(g2 x) (g2^2 + (En - u)^2))/(En - u);
B2[x_]=-((b1 E^(g1 x) g1 (g1^2 + (En - u)^2) + 
  b3 E^(g2 x) g2 (g2^2 + (En - u)^2))/((En - u) (En + u)));   

```

The constants `b1` nad `b3` to be determined using `BCs` and `g1` , `g2` are:

```
g1=Sqrt[-En^2 - u^2 + Sqrt[En^2 - u^2 + 4 En^2 u^2]];
g2=Sqrt[-En^2 - u^2 - Sqrt[En^2 - u^2 + 4 En^2 u^2]];    

```

for $x>0$ the equations become

```
-B'[x] +u A[x]==En A[x];
A'[x] +u B[x]==En B[x];    

```

solutions are

```
A[x_]=(E^(-r x) (a2 - a1 E^(2 r x)) r)/(En - u);
B[x_]=a2 E^(-r x) + a1 E^(r x);    

```

a1 and a2 constant to be determined `r=Sqrt[-(En - u)^2]`.

Finally, applying the BCS with `u = 0.125;` gives

```
MM2[W_] = 1. {Coefficient[A[0] - A1[0], {a1, a2 , b1, b3}] // Simplify,
    Coefficient[B[0] - B1[0], {a1, a2 , b1, b3}] // Simplify,
    Coefficient[B[W], {a1, a2 , b1, b3}] // Simplify,
    Coefficient[B2[0], {a1, a2 , b1, b3}] // Simplify};   

```

then we can get the eigenvalue

```
In[18]:= Eng = FindRoot[Det[MM2[20]], {En, 0}]

Out[18]= {En -> -0.0184422 + 0. I}   

```

Solving for `a1,a2,b1,b3` we get

```
{a1, a2, b1, b3}={0.404367 + 0.519605 I, -0.615275 - 0.234391 I, 
 0.045056 + 0.253905 I, -0.255964 + 0.0313093 I};   

```

this gives

```
Block[{x1 = -20, x2 = 20, En = Re@Eng[[1, 2]]}, 
 Show[Plot[{Abs[A[x]^2], (Abs@B[x])^2}, {x, 0, x2}, Frame -> True, 
       PlotRange -> {{x1, x2}, All}], 
  Plot[{(Abs@A1[x])^2, (Abs@B1[x])^2, (Abs@A2[x])^2, (Abs@
      B2[x])^2}, {x, x1, 0},  
   Frame -> True, PlotRange -> All, 
   PlotLegends -> {"A1", "B1", "A2", "B2"}]]]    

```

[![enter image description here](https://i.stack.imgur.com/gJr3C.png)](https://i.stack.imgur.com/gJr3C.png)