Modelling transparent boundary conditions on a quantum graph with Mathematica NDSolve and finite-element method

Question

Continue to this issue

I am now trying to apply a solution from the question above to three-bonded star graph. My idea is to take the third bond at [10,20] because I am using only one function u[t,x]. Link to pdetoode function:

System of equations: Method:

But on doing calculations in Mathematica, I keep getting errors and miscalculations(as for example, dimensions of the result of NDSolve are obviously off, it should correspond to number of points):

(*
Its pdetoode function, link in the post above, you can skip itenter image description here
*)
Clear[fdd, pdetoode, tooderule, pdetoae, rebuild]
fdd[{}, grid_, value_, order_, periodic_] := value;
fdd[a__] := NDSolve`FiniteDifferenceDerivative@a;

pdetoode[funcvalue_List, rest__] := 
  pdetoode[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[[1]], 
   rest];
pdetoode[{func__}[var__], rest__] := 
  pdetoode[Alternatives[func][var], rest];
pdetoode[front__, grid_?VectorQ, o_Integer, periodic_: False] := 
  pdetoode[front, {grid}, o, periodic];

pdetoode[func_[var__], time_, {grid : {__} ..}, o_Integer, 
   periodic : True | False | {(True | False) ..} : False] := 
  With[{pos = Position[{var}, time][[1, 1]]}, 
   With[{bound = #[[{1, -1}]] & /@ {grid}, 
     pat = Repeated[_, {pos - 1}], 
     spacevar = Alternatives @@ Delete[{var}, pos]}, 
    With[{coordtoindex = 
       Function[coord, 
        MapThread[
         Piecewise[{{1, PossibleZeroQ[# - #2[[1]]]}, {-1, 
             PossibleZeroQ[# - #2[[-1]]]}}, All] &, {coord, bound}]]},
      tooderule@
      Flatten@{((u : func) | 
            Derivative[dx1 : pat, dt_, dx2___][(u : func)])[x1 : pat, 
          t_, x2___] :> (Sow@coordtoindex@{x1, x2};

          fdd[{dx1, dx2}, {grid}, 
           Outer[Derivative[dt][u@##]@t &, grid], 
           "DifferenceOrder" -> o, 
           PeriodicInterpolation -> periodic]), 
        inde : spacevar :> 
         With[{i = Position[spacevar, inde][[1, 1]]}, 
          Outer[Slot@i &, grid]]}]]];

tooderule[rule_][pde_List] := tooderule[rule] /@ pde;
tooderule[rule_]@Equal[a_, b_] := 
  Equal[tooderule[rule][a - b], 0] //. 
   eqn : HoldPattern@Equal[_, _] :> Thread@eqn;
tooderule[rule_][expr_] := #[[Sequence @@ #2[[1, 1]]]] & @@ 
  Reap[expr /. rule]

pdetoae[funcvalue_List, rest__] := 
  pdetoae[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[[1]], rest];
pdetoae[{func__}[var__], rest__] := 
  pdetoae[Alternatives[func][var], rest];

pdetoae[func_[var__], rest__] := 
 Module[{t}, 
  Function[pde, #[
       pde /. {Derivative[d__][u : func][inde__] :> 
          Derivative[d, 0][u][inde, t], (u : func)[inde__] :> 
          u[inde, t]}] /. (u : func)[i__][t] :> u[i]] &@
   pdetoode[func[var, t], t, rest]]

rebuild[funcarray_, grid_?VectorQ, timeposition_: 1] := 
 rebuild[funcarray, {grid}, timeposition]

rebuild[funcarray_, grid_, timeposition_?Negative] := 
 rebuild[funcarray, grid, Range[Length@grid + 1][[timeposition]]]

rebuild[funcarray_, grid_, timeposition_: 1] /; 
  Dimensions@funcarray === Length /@ grid := 
 With[{depth = Length@grid}, 
  ListInterpolation[
     Transpose[
      Map[Developer`ToPackedArray@#["ValuesOnGrid"] &, #, {depth}], 
      Append[Delete[Range[depth + 1], timeposition], timeposition]], 
     Insert[grid, Flatten[#][[1]]["Coordinates"][[1]], 
      timeposition]] &@funcarray]

(*
constants
*)
{lb = -10, mb = 0, rb = 10, rb2 = 20, tmax = 67.3};
func2[x_] = Sin[Pi (x + 10)/10]^2

With[{u = u[t, x]}, eq = I D[u, t] + 1/2 D[u, {x, 2}] == 0;
  ic = {u == func2[x], u == 0, u == 0} /. t -> 0;
  {bcl, bcm, bcm2, bcr, 
    bcr2} = {u == 0 /. 
     x -> lb, -3 I/2 D[u, x] + D[u, t, x] + 3 I D[u, t] /. 
     x -> mb, -3 I/2 D[u, x] + D[u, t, x] + 3 I D[u, t] /. x -> rb, 
    u == 0 /. x -> rb, u == 0 /. x -> rb2}];
(*Creating grids, each corresponds to an edge of the graph
*)
points = 25; {gridl, gridr, gridr2} = 
 Array[# &, points, #] & /@ {{lb, mb}, {mb, rb}, {rb, rb2}};
difforder = 2;
(*Creating ode for each edge 
*)
{ptoofuncl, ptoofuncr, ptoofuncr2} = 
  pdetoode[u[t, x], t, #, difforder] & /@ {gridl, gridr, gridr2};

del = #[[2 ;; -2]] &;
(*Calculating ode's on both grids at each individual point 
*)
{odel, oder, oder2} = 
  del@#@eq & /@ {ptoofuncl, ptoofuncr, ptoofuncr2};
(*Calculating initial conditions on grids
*)
{odeicl, odeicr, odeicr2} = 
  MapThread[#@#2 &, {{ptoofuncl, ptoofuncr, ptoofuncr2}, ic}];
(*Calculating boundary conditions on grids
*)
{odebcl, odebcr, odebcr2} = 
  MapThread[#@#2 &, {{ptoofuncl, ptoofuncr, ptoofuncr2}, {bcl, bcr, 
     bcr2}}];
(*Calculating boundary conditions at middle point
*)
odebcm = {ptoofuncl[bcm] == ptoofuncr[bcm], 
   ptoofuncl[bcm] == ptoofuncr2[bcm2]};

odebc = {odebcm, 
   With[{sf = 1}, 
    Map[sf # + D[#, t] &, {odebcl, odebcr, odebcr2}, {2}]]};

sollst = NDSolveValue[{odel, odeicl, oder, Rest@odeicr, oder2, 
     Rest@odeicr2, odebc}, {u /@ gridl, u /@ gridr, u /@ gridr2}, {t, 
     0, tmax}, MaxSteps -> Infinity]; // AbsoluteTiming

{soll, solr, solr2} = 
  MapThread[rebuild, {sollst, {gridl, gridr, gridr2}}];


sol = {t, x} \[Function] Piecewise[{{soll[t, x], x < mb}}, solr[t, x]];
Manipulate[
 Plot[Abs[sol[t, x]]^2, {x, lb, rb}, 
  AxesLabel -> {x, 
    "|\[Psi]\!\(\*SuperscriptBox[\(|\), \(2\)]\)"}], {{t, 0, "time"}, 
  0, tmax, Appearance -> "Labeled"}]

DensityPlot[sol[t, x] // #, {t, 0, tmax}, {x, lb, rb}, 
    PlotPoints -> 50, Exclusions -> None] & /@ {Re, Im} // GraphicsRow

So is my approach correct, or I should calculate it differently?