# Solve coupled PDEs: boundary conditions problem

I need to solve 4 coupled PDEs. Thanks to @xzczd's help, a few days ago I can obtain the solution of the first 2 equations with some proper boundary conditions, see this post. Now, these two equations have been coupled with the other two, in which I encountered some new issues.

L = 30; tmax = 20; c = 1/10; \[Epsilon] = 1/10;

eqs = {D[Ar[x, t], t] + 1/4 (Ar[x, t]^3 + Ai[x, t]^2 Ar[x, t]) - D[Ar[x, t], {x, 2}] - 2 Ar[x, t] == 0,
D[Ai[x, t], t] + 1/4 (Ai[x, t]^3 + Ar[x, t]^2 Ai[x, t]) - D[Ai[x, t], {x, 2}] - 2 Ai[x, t] == 0,
D[Br[x, t], t] - D[Br[x, t], {x, 2}] + (Ar[x, t]^2 - Ai[x, t]^2)*(1/4*Br[x, t] + 10*x*D[Ar[x, t], x]
+ (50*x^2 - 1/4)* D[Ai[x, t], x]) + (Ar[x, t]^2 + Ai[x, t]^2)*(1/2*Br[x, t]
+ 10*x*D[Ar[x, t], x] - (50*x^2 + 3/4)*D[Ai[x, t], x])
+ (1/2*Bi[x, t] + 20*x*D[Ai[x, t], x] - (100*x^2 - 1/2)*D[Ar[x, t], x])*Ar[x, t]*Ai[x, t]
- 2*Br[x, t] - 2*D[Ai[x, t], {x, 3}] - 6*D[Ai[x, t], x] - 2*Ar[x, t] == 0,
D[Bi[x, t], t] - D[Bi[x, t], {x, 2}] + (Ar[x, t]^2 - Ai[x, t]^2)*(-(1/4)*Bi[x, t] - 10*x*D[Ai[x, t], x]
+ (50*x^2 - 1/4)*D[Ar[x, t], x]) + (Ar[x, t]^2 + Ai[x, t]^2)*(1/2*Bi[x, t]
+ 10*x*D[Ai[x, t], x] + (50*x^2 + 3/4)*D[Ar[x, t], x])
+ (1/2*Br[x, t] + 20*x*D[Ar[x, t], x] + (100*x^2 - 1/2)*D[Ai[x, t], x])*Ar[x, t]*Ai[x, t]
- 2*Bi[x, t] + 2*D[Ar[x, t], {x, 3}] + 6*D[Ar[x, t], x] - 2*Ai[x, t] == 0};


subject to the boundary conditions:

bc = {Ar[-L, t] == Ar[L, t], Ai[-L, t] == -Ai[L, t],
Br[-L, t] == Br[L, t], Bi[-L, t] == -Bi[L, t],
\[Epsilon]*Derivative[1, 0][Ar][-L, t] -
Ai[-L, t] == \[Epsilon]*Derivative[1, 0][Ar][L, t] - Ai[L, t],
\[Epsilon]*Derivative[1, 0][Ai][-L, t] +
Ar[-L, t] == -\[Epsilon]*Derivative[1, 0][Ai][L, t] - Ar[L, t],
\[Epsilon]*Derivative[1, 0][Br][-L, t] -
Bi[-L, t] == \[Epsilon]*Derivative[1, 0][Br][L, t] - Bi[L, t],
\[Epsilon]*Derivative[1, 0][Bi][-L, t] +
Br[-L, t] == -\[Epsilon]*Derivative[1, 0][Bi][L, t] - Br[L, t]};


The initial conditions are set randomly:

ic = {Ar[x, 0] == iniAr[x], Ai[x, 0] == iniAi[x], Br[x, 0] == iniBr[x], Bi[x, 0] == iniBi[x]};

iniAr = ListInterpolation[RandomReal[{-c, c}, 20], {{-L, L}}];
iniAi = ListInterpolation[RandomReal[{-c, c}, 20], {{-L, L}}];
iniBr = ListInterpolation[RandomReal[{-c, c}, 20], {{-L, L}}];
iniBi = ListInterpolation[RandomReal[{-c, c}, 20], {{-L, L}}];


As mentioned in the previous post, when solving it with NDSolve

sol = NDSolveValue[{eqs, ic, bc}, {Ar[x, t], Ai[x, t], Br[x, t], Bi[x, t]}, {x, -L, L}, {t, 0, tmax}]


it returns an error

NDSolveValue::bcedge: Boundary condition Ai[-30,t]==-Ai[30,t] is not specified on a single edge of the boundary of the computational domain.

So I tried @xzczd's pdetoode

points = 200; domain = {-L, L}; difforder = 4;
grid = Array[# &, points, domain];
ptoofunc = pdetoode[{Ar, Ai, Br, Bi}[x, t], t, grid, difforder];
odebc = Map[ptoofunc, bc, {2}]

del = #[[2 ;; -2]] &;
odeic = del /@ ptoofunc@ic;
ode = del /@ ptoofunc@eqs;
sollst = NDSolveValue[{ode, odeic, odebc}, Table[v[x], {v, {Ar, Ai, Br, Bi}}, {x, grid}], {t, 0, tmax}];


However, I got another error:

NDSolveValue::mconly: For the method NDSolveIDA, only machine real code is available. Unable to continue with complex values or beyond floating-point exceptions.

What is machine real code? Could anyone help with solving the problem and explaining the error? Thank you in advance!

Strongly related:

Why does NDSolve fail to solve the PDEs and spit out mconly warning?

odebc = With[{sf = 1}, Map[ptoofunc[sf # + D[#, t]] &, bc, {2}]];

odeic = ptoofunc@ic;

var = Table[v[x], {v, {Ar, Ai, Br, Bi}}, {x, grid}];

lhs = D[Flatten[var][t] // Through, t];

{barray, marray} =
CoefficientArrays[Flatten@{ode, odebc}, lhs]; // AbsoluteTiming
(* {0.118406, Null} *)

rhs = LinearSolve[marray, -barray]; // AbsoluteTiming
(* {2.87564, Null} *)

sollst = NDSolveValue[{lhs == rhs, odeic}, var, {t, 0, tmax}]; // AbsoluteTiming
(* {3.15959, Null} *)
`

Necessary explanation has been included in the post above so I'd like not to repeat it here. (Don't miss those links in that post. )