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?