# Solving PDE for a parametric function, using ParametricNDSolve

I've to solve a $4 \times 4$ matrix PDE which involves the following matrices:

Mp[r_] = {{-kx Sqrt[1/r^2] - Sqrt[r^4]/r^2 + Sqrt[r^2] (1 - 1/r + ω),
-2 I (1 - 1/r), -ky Sqrt[1/r^2], 0}, {2 I (1 - 1/r), kx Sqrt[1/r^2]
+ Sqrt[r^4]/r^2 - Sqrt[r^2] (1 - 1/r - ω), 0, -ky Sqrt[1/r^2]},
{ky Sqrt[1/r^2], 0, -kx Sqrt[1/r^2] - Sqrt[r^4]/r^2 -
Sqrt[r^2] (1 - 1/r + ω), -2 I (1 - 1/r)},
{0, ky Sqrt[1/r^2], 2 I (1 - 1/r), kx Sqrt[1/r^2] +
Sqrt[r^4]/r^2 + Sqrt[r^2] (1 - 1/r - ω)}}

Mn[r_] = {{-kx Sqrt[1/r^2] - Sqrt[r^4]/r^2 - Sqrt[r^2] (1 - 1/r + ω),
-2 I (1 - 1/r), ky Sqrt[1/r^2], 0}, {2 I (1 - 1/r), kx Sqrt[1/r^2] +
Sqrt[r^4]/r^2 + Sqrt[r^2] (1 - 1/r - ω), 0, ky Sqrt[1/r^2]},
{-ky Sqrt[1/r^2], 0, -kx Sqrt[1/r^2] - Sqrt[r^4]/r^2 + Sqrt[r^2] (1
- 1/r + ω), -2 I (1 - 1/r)}, {0, -ky Sqrt[1/r^2], 2 I (1 -1/r),
kx Sqrt[1/r^2] + Sqrt[r^4]/r^2 - Sqrt[r^2] (1 - 1/r - ω)}}

kk = √(kx^2 + ky^2);

GIR = DiagonalMatrix[{ -( √(kk + ω)/√(kk - \
ω)) ,  √(kk - ω)/√(kk + ω),  \
√(kk - ω)/√(kk + ω), -( √(kk + \
ω)/√(kk - ω)) }]


I'm trying solve for a 3-parameter (kx, ky, ω) function $G$ using the following PDE. Note GIR defined above gives a boundary condition for $r$ close to zero, since $r=0$ is singularity I've kept it $0.001$). The PDE is

soln = ParametricNDSolve[{G'[r] == Mp[r] - G[r].Mn[r].G[r], G[0.001] == GIR},
G, {r, 1, 10}, {kx, ky, ω}]
G[kx, ky, ω][r] /. soln


After running this code, when I evaluate e.g. G[2, 1, 1, 0.001] I get back G[2, 1, 1, 0.001] ; I can't figure out why $G$ is not being evaluated. In fact the final quantity of my interest is:

M[r_, kx_, ky_, ω_] := Tr[Im[Evaluate[G[kx, ky, ω][r] /. soln]]]
ContourPlot[M[r, kx, 0, 0], {kx, 0, 5}, {r, 1, 10}, FrameLabel -> {"k", "Tr[Im M]"}]


But this gives an empty plot. Any kind of help will be appreciated, thanks!

G[r_] := Array[g[##][r] &, {4, 4}]
soln = ParametricNDSolve[
Join[Thread[ Flatten@Array[g[##]'[r] &, {4, 4}] == Flatten[Mp[r] -G[r].Mn[r].G[r]]],
Flatten[Array[g[##] &, {4, 4}]],
{r, .001, 10}, {kx, ky, ω}]

g[1, 2][1, 2, 3][1] /. soln
(*-3.58322*10^-17 + 0.22098 I*)


The following is probably numerical noise:

M[r_?NumericQ, kx_?NumericQ, ky_?NumericQ, ω_?NumericQ] :=
Tr[Re[Evaluate[Array[g[##][kx, ky, ω][r] &, {4, 4}] /. soln]]]
ContourPlot[M[r, kx, 0, 0], {kx, 0, 5}, {r, 1, 10},
FrameLabel -> {"k", "Tr[Im M]"}]


• Thanks, that was really helpful yet I still get a blank graph only. Nov 10, 2015 at 2:09
• @Skylar15 There are other problems. For certain parameters your system is stiff. For others the solution involves complex numbers.My code is just a start for you to be able to get some numerical output Nov 10, 2015 at 2:11
• @Skylar15 See edit Nov 10, 2015 at 2:25
• Thanks! Yeah, I am trying to take care of the singularities and numerical noises. For those I've to fix some other parameters. Appreciate your effort! Nov 10, 2015 at 2:36