# Error using pdetoode

I am trying to solve the system of equations:

eq1 = (Exp[-2 λ[t, r]] D[R[t, r], r]^2)/R[t, r]^2 -
2 (Exp[-2 λ[t, r]] D[R[t, r], r] D[λ[t, r], r])/
R[t, r] +
2 (Exp[-2 λ[t, r]] D[R[t, r], {r, 2}])/R[t, r] - (
Exp[-2 ν[t, r]] D[R[t, r], t]^2)/R[t, r] -
2 (Exp[-2 ν[t, r]] D[R[t, r], t] D[ν[t, r], t])/R[t, r] + (
8 Pi G)/3 ϵ[t, r];

eq3 = (Exp[-2 λ[t, r]] D[R[t, r], r]^2)/R[t, r]^2 +
2 (Exp[-2 λ[t, r]] D[R[t, r], r] D[ν[t, r], r])/
R[t, r] - (Exp[-2 ν[t, r]] D[R[t, r], t]^2)/R[t, r]^2 +
2 (Exp[-2 ν[t, r]] D[R[t, r], t] D[ν[t, r], t])/R[t, r] -
2 (Exp[-2 ν[t, r]] D[R[t, r], {t, 2}])/R[t, r] - (8 Pi G)/
3 ϵ[t, r]/3;


where the functions $$\epsilon$$ and $$\nu$$ are given by:

ϵ[t_, r_]=(E^( - 4/3 λ[t, r]) (1/100)^(8/3))/R[t, r]^(8/3)
ν[t_, r_]=1/4 (Log[1/(1 + t)^(8/3)] - Log[E^(-(4/3) λ[t, r])/R[t, r]^(8/3)])


I use the method pdetoode by @xzczd as follows:

lb = 10^-1;
rb = 1;
xdifforder1 = 2;
xdifforder2 = 1;
ic1 = {R[0, r] == 10^-2, Derivative[1, 0][R][0, r] == 1};
ic2 = {λ[0, r] == 0};
bc1 = {R[t, 0] == t + 10^-2, R[t, 10^-1] == 10^-2,
Derivative[0, 1][R][t, 0] == 0};
bc2 = {λ[t, 0] == 0};
points = 25;
grid = Array[# &, points, {lb, rb}];
removeredundant1 = #[[2 ;; -2]] &;
removeredundant2 = #[[1 ;; -1]] &;
ptoofunc1 = pdetoode[R[t, r], t, grid, xdifforder1];
ptoofunc2 = pdetoode[λ[t, r], t, grid, xdifforder2];
odeqn1 = eq1 // ptoofunc1 // removeredundant1;
odeqn2 = eq3 // ptoofunc2 // removeredundant2;
odeic1 = removeredundant1 /@ ptoofunc1@ic1;
odeic2 = removeredundant2 /@ ptoofunc2@ic2;
odebc1 = bc1 // ptoofunc1;
odebc2 = bc2 // ptoofunc2;
sollst = NDSolveValue[{odebc1, odebc2, odeic1, odeic2, odeqn1, odeqn2}, {R /@ grid, [Lambda] /@ grid}, {t, 0, 1},  MaxSteps -> Infinity]
sol = rebuild[sollst, grid]


but I get the error:

So it seems that while it works well for $$R$$, the grid points are not mapped onto $$\lambda$$ and it becomes a function of a list. I do not understand why it happens, since it seems to me that I am applying the same method to both $$R$$ and $$\lambda$$. Can anyone help?

EDIT

After the comment by @xzczd, my code now looks like this:

lb = 10^-1;
rb = 1;
xdifforder = 4;
ic = {R[0, r] == 10^-2,
Derivative[1, 0][R][0, r] == 1, λ[0, r] == 0};
bc = {R[t, 0] == t + 10^-2, R[t, 10^-1] == 10^-2,
Derivative[0, 1][R][t, 0] == 0, λ[t, 0] == 0};
points = 25;
grid = Array[# &, points, {lb, rb}];
removeredundant = #[[3 ;; -3]] &;
ptoofunc = pdetoode[{R[t, r], λ[t, r]}, t, grid, xdifforder];
odeqn1 = eq1 // ptoofunc // removeredundant;
odeqn2 = eq3 // ptoofunc // removeredundant;
odeic = removeredundant /@ ptoofunc@ic;
odebc = bc // ptoofunc;


but now Mathematica complaints that there are fewer variables than equations

• Why do you use pdetoode twice? Commented May 20, 2021 at 7:02
• @xzczd because I have two functions $R$ and $\lambda$, I have found some posts on this site for system of equations and, as far as I understand, this seems to be the method. Is it wrong? Commented May 20, 2021 at 7:05
• Yes, it's wrong. If the domain of the unknown functions are the same, then you only need to use pdetoode once, i.e. pdetoode[{R[t, r], \[Lambda][t, r]}, … or simply pdetoode[{R, \[Lambda]}[t, r], … Commented May 20, 2021 at 7:10
• @xzczd but the functions have different differentiation order: $R$ it is 2, while $\lambda$ it is 1 Commented May 20, 2021 at 7:13
• diff is short for difference. You may want to read the comments here: mathematica.stackexchange.com/a/174775/1871 Commented May 20, 2021 at 7:21

This system can be solved with DAE solver on the short time interval {t,0,.08} as follows

ϵ = (Exp[(-4/3 λ[t, r])] (1/100)^(8/3))/R[t, r]^(8/3);
ν = 1/
4 (Log[1/(1 + t)^(8/3)] -
Log[E^(-(4/3) λ[t, r])/R[t, r]^(8/3)]);

eq = {(Exp[-2 λ[t, r]] D[R[t, r], r]^2)/R[t, r]^2 -
2 (Exp[-2 λ[t, r]] D[R[t, r], r] D[λ[t, r], r])/
R[t, r] +
2 (Exp[-2 λ[t, r]] D[R[t, r], {r, 2}])/
R[t, r] - (Exp[-2 ν] D[R[t, r], t]^2)/R[t, r] -
2 (Exp[-2 ν] D[R[t, r], t] D[ν, t])/R[t, r] + (8 Pi G)/
3 ϵ,
(Exp[-2 λ[t, r]] D[R[t, r], r]^2)/R[t, r]^2 +
2 (Exp[-2 λ[t, r]] D[R[t, r], r] D[ν, r])/
R[t, r] - (Exp[-2 ν] D[R[t, r], t]^2)/R[t, r]^2 +
2 (Exp[-2 ν] D[R[t, r], t] D[ν, t])/R[t, r] -
2 (Exp[-2 ν] D[R[t, r], {t, 2}])/R[t, r] - (8 Pi G)/
3 ϵ/3} // Simplify

ic = {R[t, r] == 10^-2, D[R[t, r], t] == 1, λ[t, r] == 0} /.
t -> 0;
bc = {R[t, r] - t - 10^-2 == 0,
D[R[t, r], r] == 0, λ[t, r] == 0} /. r -> lb;

G = 1; lb = 10^-1;
rb = 1;
sol = NDSolve[{eq[[1]]+eq[[2]] == 0, eq[[2]] == 0, ic,
bc}, {R, λ}, {t, 0, .08}, {r, lb, rb},
Method -> {"IndexReduction" -> Automatic,
"EquationSimplification" -> "Residual",
"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 80, "MaxPoints" -> 80,
"DifferenceOrder" -> 8}}}]

{Plot3D[Evaluate[R[t, r] /. sol[[1]]], {t, 0, 0.08}, {r, lb, rb},
ColorFunction -> Hue, Mesh -> None, AxesLabel -> Automatic,
PlotRange -> All, Boxed -> False, PlotPoints -> 50],
Plot3D[Evaluate[λ[t, r] /. sol[[1]]], {t, 0, 0.08}, {r, lb,
rb}, ColorFunction -> Hue, Mesh -> None, AxesLabel -> Automatic,
Boxed -> False, PlotPoints -> 50]}


It looks very regular, but it is not clear why system blow up after t=0.089

There is a stable solution for $$\lambda >0$$ with the set of boundary conditions

ic = {R[t, r] == 10^-2,
D[R[t, r], t] == 1/100, λ[t, r] == 0} /. t -> 0;
bc = {R[t, r] - t/100 - 10^-2 == 0,
D[R[t, r], r] == 0, λ[t, r] == 0} /. r -> lb;


This solution on {t,0,1} looks like

Note, we are not been able to reproduce any of these results with using pdetoode or pdetoae.