I often solve pdes for my research, and years ago I found pdetoode in this forum is very handy. Although it is a small piece of code, it solves several interesting and challenging pdes (here, here and here for example) in this forum.

Thus, I decied to figure out its working principle, because I strongly believe that it is greatly benificial to both ode and Mathematica understandings.

I have been studying hard on the package pdetoode by xzczd for many times and many days each time I picks it up.

I still can't figure out the how the pde converted to odes.

In particularlly, the following command:

((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])

I guess from the context that pat here means whatever repeats itself several times exactly. However, for code here, after I print dx1, and x1. In some cases dx1 = 0 1, x1=x y , this confues me a lot.

I also attempt to use Trace on pdetoode, but its long long out confuses me more.

How can this be? Maybe I misunderstand something here? What does this piece of code try to do?

I tried my best but I still fails to understand it.

Can anyone explain something on the above code or pdetoode? Thanks!

In order to analyze how pdetoode works, below is the code that I used to Print some local variables. The confusing printed information mentioned above is shown in the following figure.

Remove["Global`*"] // Quiet;

Clear[fdd, pde2ode, tooderule, pdetoae, rebuild]

(*====== fdd =======*)

fdd[{}, grid_, value_, order_, periodic_] := value;
fdd[a__] := NDSolve`FiniteDifferenceDerivative@a;


(*====== pde2ode =======*)

pde2ode[funcvalue_List, rest__] := 
  With[{s = Style[#, Purple, Bold] &}, 
   Print[s@"pde2ode-1 ********************************"];
   Print[s@"funcvalue = ", s@funcvalue]; Print[s@"rest = ", s@rest];
   pde2ode[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[[1]], 
    rest]];

pde2ode[{func__}[var__], rest__] := 
  With[{s = Style[#, Red, Bold] &}, 
   Print[s@"pde2ode-2 ********************************"];
   Print[s@"func = ", s@func]; Print[s@"var = ", s@var]; 
   Print[s@"rest = ", s@rest];
   pde2ode[Alternatives[func][var], rest]];

pde2ode[front__ (* front : u[x,t],t *), grid_?VectorQ, o_Integer, 
   periodic_ : False] := 
  With[{s = Style[#, Blue, Bold] &}, 
   Print[s@"pde2ode-3 ********************************"];
   Print[s@"front = ", s@front]; Print[s@"grid = ", s@grid];
   pde2ode[front, {grid}, o, periodic]];

pde2ode[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}]]},
     With[{s = Style[#, Darker[Green], Bold] &}, 
      Print[s@"pde2ode-core *****************************"];
      Print[s@"func = ", s@func]; Print[s@"time = ", s@time]; 
      Print[s@"grid = ", s@grid];
      Print[s@"pos = ", s@pos]; Print[s@"bound = ", s@bound]; 
      Print[s@"spacevar = ", s@spacevar]; 
      Print[s@"coordtoindex = ", s@coordtoindex]];
     tooderule@Flatten@{
        (*------- 
        rule_1 --------*)
        ((u : func) | 
            Derivative[dx1 : pat, dt_, dx2___][(u : func)])[x1 : pat, 
          t_, 
          x2___] :>
         (With[{}, Print["-------------"]; 
           Print["u = ", u]; Print["dx1 = ", dx1]; 
           Print["pat = ", pat]; Print["dt = ", dt]; 
           Print["dx2 = ", dx2]; Print["x1 = ", x1]; 
           Print["x2 = ", x2]];
          Sow@coordtoindex@{x1, x2}; 
          With[{}, Print["coord = ", {x1, x2}]; 
           Print["coordtoindex[coord] = ", coordtoindex@{x1, x2}]];
          
          fdd[{dx1, dx2}, {grid}, 
           Outer[Derivative[dt][u@##]@t &, grid], 
           "DifferenceOrder" -> o, PeriodicInterpolation -> periodic]),
        (*------- rule_2 --------*)
        
        inde : spacevar :> 
         With[{i = Position[spacevar, inde][[1, 1]]}, 
          Outer[Slot@i &, grid]]
        }]]];


(*====== tooderule =======*)

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 =======*)

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]] &@
   pde2ode[func[var, t], t, rest]]


(*====== rebuild =======*)

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]
(*=================== example of pdetoode ================*)

L = 1; T = 1; x0 = -L/4; sigma = L/30;
domain = {-L/2, L/2};
{eq1, eq2} = {D[G[x, y, t], t] == -(D[G[x, y, t], x] + D[G[x, y, t], y]) - 
     I (f[x, y, t] + f2[x, y, t]), 
   D[f[x, t], t] == -D[f[x, t], x] + f[x, t] - I (G2[x, t])};

{ic1, ic2} = {G[x, y, 0] == 
    Exp[-((x - x0)/(Sqrt[2] sigma))^2 - ((y - x0)/(Sqrt[2] sigma))^2], 
   f[x, 0] == 0};

{bc1, bc2} = {G[x, y, t] == 0 /. Outer[{# -> #2} &, {x, y}, domain], 
   f[x, t] == 0 /. List /@ Thread[x -> domain]};

points = 5;
grid = Array[# &, points, domain];

difforder = 2;

ptoofunc1 = 
 pde2ode[{G[x, y, t], f[x, y, t], f2[x, y, t]}, t, {grid, grid}, difforder]

grid
ptoofunc2 = pde2ode[{G2[x, t], f[x, t]}, t, grid, difforder];

ptoofunc1 = 
 pde2ode[{G, f, f2}[x, y, t], t, {grid, grid}, difforder]

ptoofunc2 = pde2ode[{G2, f}[x, t], t, grid, difforder];


del = #[[2 ;; -2]] &;
rule1 = {f[x_, y_][t_] :> f[x][t], f2[x_, y_][t_] :> f[y][t]};
ode1 = del /@ del@ptoofunc1@eq1;
ode1 = del /@ del@ptoofunc1@eq1 /. rule1;

rule2 = G2[x_][t_] :> G[x, 0][t];
ode2 = del@ptoofunc2@eq2 /. rule2;

odeic1 = ptoofunc1@ic1;

odeic2 = ptoofunc2@ic2;


diff = With[{sf = 1}, D[#, t] + #] &;
odebc1 = Map[diff, MapAt[del /@ # &, ptoofunc1@bc1, {1}], {-2}];

odebc2 = Map[diff, ptoofunc2@bc2, {-2}];

sol = NDSolveValue[{ode1, ode2, odeic1, odeic2, odebc1, 
    odebc2}, {Outer[G, grid, grid], f /@ grid}, {t, 0, T}];

solG = rebuild[sol[[1]], {grid, grid}, 3];
solf = rebuild[sol[[2]], grid, 2];

(*Manipulate[Plot3D[solG[x,y,t]//Evaluate,{x,##},{y,##},PlotRange->{-0.1,1}],{\
t,0,T}]&@@domain*)

Manipulate[
   Plot[solf[x, t] // Abs // Evaluate, {x, ##}, PlotRange -> {0, 0.2}], {t, 0,
     T}] & @@ domain