I am trying to solve a 1D non-linear ODE with FEM. The ODE corresponds to a Model with 3 parameters ([Beta], [lambda], A_2). The ODE reads
mpl/(m^2)Cos(z)^4 f''(z)-2mpl/(m^2)Cos(z)^4Tan(z)f'(z)-A2\rho(z)/mpl *f(z) + Exp(-[lambda]f(z)/mpl+Log[[lambda]+Log[10][Beta]),
where mpl, m are some fixed constants, \rho(z) is a density. I have implemented the above equation and Mathematica does return a solution for the parameters that are of interest to me. I give an initial seed to NDSolve, and create my own mesh.
For many parameters it does return the correct solution (I have a very good analytical approximation to check that) but for some it hardly improves on the not very accurate seed and still returns a solution.
My question is: For the parameters [beta]=1, A2=10^(40), \lambda=10^27; it does return the correct solution, for [beta]=1, A2=10^(40), \lambda=10^(31) it does not and seems to get stuck on the seed. What is the reason for this? I provide the full code necessary to just run the final differential equation with NDSolveValue. The following cell has to be activated first (constants definition, seed definition, mesh definition, ODE definition etc):
Needs["NDSolve`FEM`"];
Needs["NumericalCalculus`"]
(* First define conversion factors between various units*)
CPrec = 500; (* numerical precission for later calculations*)
constants =
SetPrecision[{hbar -> 1.05457*^-34, c0 -> 2.997925*^8,
ec -> 1.602177*^-19}, 5 CPrec];
(* Note that ec is the unitless(!!) value of the electric charge, \
which is used to convert eV to J. ev=ec*J *)
RES = SetPrecision[149597870700*m2invMeV,
5 CPrec]; (*Earth-Sun distance also known as AE or AU, \
Subscript[r, AU] in the paper*)
aG = SetPrecision[1.30199*10^(-32),
5 CPrec]; (*acceleration of earth towards sun in MeV, Pitschmann \
value*)
GN = SetPrecision[6.70861*^-45,
5 CPrec]; (*newtons constant in MeV^(-2), Pitschmann value*)
mpl = SetPrecision[2.4353635268*^21, 5 CPrec];(*planck mass in MeV*)
VC = SetPrecision[6.73304*^-34,
5 CPrec]; (*in MeV^4, this is the effective potential for the \
vacuum density and the vaccuum expactation value of \[Phi]*)
m2invMeV =
1/(hbar c0/(1*^6*ec)) /.
constants; (*to be read as convert meter to MeV^-1*)
kg2MeV = 1/(ec 10^6/(c0^2)) /. constants;
s2invMeV = 1/(hbar/(1*^6 ec)) /. constants;
N2invMeV2 = kg2MeV m2invMeV/s2invMeV^2; (*Newton to inverse MeV^2*)
kgm32MeV4 =
kg2MeV/m2invMeV^3; (*to be read as kg/m^3 converted to MeV^4 / \
density to natural units*)
(*Next define various constants / parameters that are use in the LLR \
analysis*)
rhoV = SetPrecision[1*^-15, 5 CPrec]; (* vacuum density in MeV^4 *)
rhoE = SetPrecision[0.0000236851,
5 CPrec]; (*mean density of the earth in MeV^4 *)
rhoS = SetPrecision[1408*kgm32MeV4,
5 CPrec]; (* mean density of the sun in MeV^4 *)
rhoM = SetPrecision[0.0000143877,
5 CPrec]; (*mean density of the moon in MeV^4, Pitschmannvalue*)
rhoN = SetPrecision[1.435550104*10^10,
5 CPrec]; (*neutron density in MeV^4, Pitschmannvalue*)
RE = SetPrecision[6378100*m2invMeV,
5 CPrec]; (*equatorial radius of the earth in MeV*)
RS = SetPrecision[695700000*m2invMeV,
5 CPrec]; (*equatorial radius of the sun in MeV*)
RM = SetPrecision[1738100*m2invMeV,
5 CPrec]; (*equatorial radius of the moon in MeV*)
RN = SetPrecision[0.0025,
5 CPrec]; (*Neutron radius MeV^-1 - Pitschmann value*)
ME = SetPrecision[5.9724*10^(24)*kg2MeV,
5 CPrec]; (*total mass of the earth in MeV*)
MS = SetPrecision[1.9891*10^(30)*kg2MeV,
5 CPrec]; (*total mass of the sun in MeV*)
MN = SetPrecision[939.565346,
5 CPrec]; (* neutron mass in MeV, Pitschmann value *)
REM = SetPrecision[385000558.4*m2invMeV,
5 CPrec]; (*Earth-Moon distance in MeV^-1*)
REMax = SetPrecision[405000000*m2invMeV,
5 CPrec]; (*Max. Earth-Moon distance in MeV^-1*)
\[Rho]V = rhoV;
\[Rho]E = rhoE;
\[Rho]S = rhoS;
\[Rho]M = rhoM;
\[Rho]N = rhoN;
\[Rho]M2 =
SetPrecision[2514.0`*kgm32MeV4, 5 CPrec]; (*mirror density cannex*)
ddr = SetPrecision[5.0677398985414185`*^7, 1000];
phirho[A_, \[Lambda]_, \[Beta]_, \[Rho]_] :=
SetPrecision[
mpl/\[Lambda] If[
2*Log10[\[Lambda]] + \[Beta] - Log10[A] - Log10[\[Rho]] < 1000000,
ProductLog[(\[Lambda]^2 10^\[Beta])/(A \[Rho])],
Log[\[Lambda]^2/(A \[Rho])] +
1/(6 (\[Beta] Log[10] +
Log[\[Lambda]^2/(A \[Rho])])^3) (6 \[Beta] Log[
10] (\[Beta] Log[10] + Log[\[Lambda]^2/(A \[Rho])])^3 -
6 (-1 + \[Beta] Log[10] +
Log[\[Lambda]^2/(A \[Rho])]) (1 + \[Beta]^2 Log[10]^2 +
2 \[Beta] Log[10] Log[\[Lambda]^2/(A \[Rho])] +
Log[\[Lambda]^2/(A \[Rho])]^2) Log[\[Beta] Log[10] +
Log[\[Lambda]^2/(A \[Rho])]] +
3 (-3 + \[Beta] Log[10] +
Log[\[Lambda]^2/(A \[Rho])]) Log[\[Beta] Log[10] +
Log[\[Lambda]^2/(A \[Rho])]]^2 +
2 Log[\[Beta] Log[10] + Log[\[Lambda]^2/(A \[Rho])]]^3)],
1000]
murho[\[Lambda]_, \[Rho]_, A2_] :=
SetPrecision[
1/mpl Sqrt[\[Lambda]^2 V0[A2, \[Lambda]] Exp[
SetPrecision[-\[Lambda] phirho[\[Lambda], \[Rho], A2]/mpl,
5 CPrec]] + A2 \[Rho]], 5 CPrec];
Veff[\[Lambda]_, \[Rho]_, \[Phi]_, A2_] :=
SetPrecision[
V0[A2, \[Lambda]] Exp[
SetPrecision[-\[Lambda] \[Phi]/mpl,
5 CPrec]] + (A2 \[Rho])/(2 mpl^2) \[Phi]^2, 5 CPrec];
j[x_] := 3/x^2 (1 - Tanh[x]/x)
js[x_] :=
1 - (2 x^2)/5 + (17 x^4)/105 - (62 x^6)/945 + (1382 x^8)/51975
SC[\[Lambda]_, \[Rho]_, A2_, R_] :=
If[murho[\[Lambda], \[Rho], A2] R <= 0.01,
SetPrecision[
js[murho[\[Lambda], \[Rho], A2] R]/(1 +
murho[\[Lambda], rhoV, A2]/murho[\[Lambda], \[Rho], A2] Tanh[
murho[\[Lambda], \[Rho], A2] R]), 5 CPrec],
SetPrecision[
j[murho[\[Lambda], \[Rho], A2] R]/(1 +
murho[\[Lambda], rhoV, A2]/murho[\[Lambda], \[Rho], A2] Tanh[
murho[\[Lambda], \[Rho], A2] R]),
5 CPrec]] (*screening charge*)
Q[\[Lambda]_, \[Rho]_, A2_, R_] := SC[\[Lambda], \[Rho], A2, R]
\[Mu][A_, \[Lambda]_, \[Beta]_, \[Rho]_] :=
SetPrecision[
Sqrt[SetPrecision[\[Lambda]^2 E^
SetPrecision[-\[Lambda] \[Phi]M[
A, \[Lambda], \[Beta], \[Rho]]/mpl + \[Beta]/Log10[E],
1000] + A \[Rho], 1000]]/mpl, 1000]
\[Phi]M[A_, \[Lambda]_, \[Beta]_, \[Rho]_] :=
phirho[A, \[Lambda], \[Beta], \[Rho]]
\[Mu]0[\[Lambda]_?NumberQ,
A2_?NumberQ, \[Beta]_?NumberQ, \[Phi]0_?NumberQ, \[Rho]V2_?
NumberQ] :=
SetPrecision[
1/mpl Sqrt[\[Lambda]^2 E^
SetPrecision[-\[Lambda] \[Phi]0/mpl + \[Beta]/Log10[E], 1000] +
A2 \[Rho]V2], 1000]
D0[\[Lambda]_?NumberQ,
A2_?NumberQ, \[Beta]_?NumberQ, \[Phi]0_?NumberQ, \[Rho]V2_?
NumberQ] :=
SetPrecision[\[Lambda] /mpl E^
SetPrecision[-\[Lambda] \[Phi]0/mpl + \[Beta]/Log10[E],
1000] - (A2 \[Rho]V2)/mpl^2 \[Phi]0, 1000]
\[Phi]d[\[Lambda]_, A2_, \[Beta]_, \[Phi]0_, d_, \[Rho]V2_] :=
SetPrecision[\[Phi]0 +
D0[\[Lambda],
A2, \[Beta], \[Phi]0, \[Rho]V2]/\[Mu]0[\[Lambda],
A2, \[Beta], \[Phi]0, \[Rho]V2]^2 (1 -
Cosh[SetPrecision[\[Mu]0[\[Lambda],
A2, \[Beta], \[Phi]0, \[Rho]V2]*d, 5 CPrec]]), 5 CPrec]
\[Phi]2[\[Lambda]_, A2_, \[Beta]_, \[Phi]0_, d_,
z_, \[Rho]M2_, \[Rho]V2_] :=
SetPrecision[
Piecewise[{{\[Phi]0 +
D0[\[Lambda],
A2, \[Beta], \[Phi]0, \[Rho]V2]/\[Mu]0[\[Lambda],
A2, \[Beta], \[Phi]0, \[Rho]V2]^2 (1 -
Cosh[SetPrecision[\[Mu]0[\[Lambda],
A2, \[Beta], \[Phi]0, \[Rho]V2]*(z), 5 CPrec]]),
d - Abs[z] > 0},
{\[Phi]M[
A2, \[Lambda], \[Beta], \[Rho]M2] + (\[Phi]d[\[Lambda],
A2, \[Beta], \[Phi]0, d, \[Rho]V2] - \[Phi]M[
A2, \[Lambda], \[Beta], \[Rho]M2]) E^
SetPrecision[-\[Mu][
A2, \[Lambda], \[Beta], \[Rho]M2] (Abs[z] - d), 5 CPrec],
Abs[z] - d > 0}}], 5 CPrec]
\[Phi]\[Rho][A_, \[Lambda]_, \[Beta]_, \[Rho]_] :=
SetPrecision[
mpl/\[Lambda] If[
2*Log10[\[Lambda]] + \[Beta] - Log10[A] - Log10[\[Rho]] < 1000000,
ProductLog[(\[Lambda]^2 10^\[Beta])/(A \[Rho])],
Log[\[Lambda]^2/(A \[Rho])] +
1/(6 (\[Beta] Log[10] +
Log[\[Lambda]^2/(A \[Rho])])^3) (6 \[Beta] Log[
10] (\[Beta] Log[10] + Log[\[Lambda]^2/(A \[Rho])])^3 -
6 (-1 + \[Beta] Log[10] +
Log[\[Lambda]^2/(A \[Rho])]) (1 + \[Beta]^2 Log[10]^2 +
2 \[Beta] Log[10] Log[\[Lambda]^2/(A \[Rho])] +
Log[\[Lambda]^2/(A \[Rho])]^2) Log[\[Beta] Log[10] +
Log[\[Lambda]^2/(A \[Rho])]] +
3 (-3 + \[Beta] Log[10] +
Log[\[Lambda]^2/(A \[Rho])]) Log[\[Beta] Log[10] +
Log[\[Lambda]^2/(A \[Rho])]]^2 +
2 Log[\[Beta] Log[10] + Log[\[Lambda]^2/(A \[Rho])]]^3)], 1000]
EQN\[Phi]0[\[Lambda]_, A2_, \[Beta]_, \[Phi]0_, \[Rho]_, \[Rho]V2_,
d_] := SetPrecision[-\[Phi]0 + ((mpl \[Lambda] -
A2 E^SetPrecision[(\[Lambda] \[Phi]0)/mpl - \[Beta] Log[10],
1000] \[Rho]V2 \[Phi]0) (-1 +
Cosh[SetPrecision[(d Sqrt[
E^SetPrecision[-((\[Lambda] \[Phi]0)/mpl) + \[Beta] Log[
10], 1000] \[Lambda]^2 + A2 \[Rho]V2])/mpl,
1000]]))/(\[Lambda]^2 +
A2 E^SetPrecision[(\[Lambda] \[Phi]0)/mpl - \[Beta] Log[10],
1000] \[Rho]V2) + \[Phi]\[Rho][
A2, \[Lambda], \[Beta] , \[Rho]] + ((E^
SetPrecision[-((\[Lambda] \[Phi]0)/mpl) + \[Beta] Log[10],
1000] mpl \[Lambda] -
A2 \[Rho]V2 \[Phi]0) Sinh[(d Sqrt[
E^SetPrecision[-((\[Lambda] \[Phi]0)/mpl) + \[Beta] Log[
10], 1000] \[Lambda]^2 + A2 \[Rho]V2])/mpl])/(Sqrt[
E^SetPrecision[-((\[Lambda] \[Phi]0)/mpl) + \[Beta] Log[10],
1000] \[Lambda]^2 + A2 \[Rho]V2] Sqrt[
A2 \[Rho] (1 + (\[Lambda] \[Phi]\[Rho][
A2, \[Lambda], \[Beta] , \[Rho]])/mpl)]), 1000]
\[Phi]0[\[Lambda]_?NumericQ, A2_, \[Beta]_, d_, \[Rho]_, \[Rho]V2_] :=
Block[{res, startp},
startp =
Abs[SetPrecision[\[Phi]\[Rho][A2, \[Lambda], \[Beta], \[Rho]V2],
5 CPrec]];
res =
N[(phi0 /.
SetPrecision[
FindRoot[
Re[SetPrecision[
EQN\[Phi]0[SetPrecision[\[Lambda], 1000],
SetPrecision[A2, 1000], SetPrecision[\[Beta], 1000],
SetPrecision[phi0, 1000], SetPrecision[\[Rho], 1000],
SetPrecision[\[Rho]V2, 1000], SetPrecision[d, 1000]],
5 CPrec]], {phi0, startp, 0,
SetPrecision[1.5, 5 CPrec] startp}, AccuracyGoal -> 5 CPrec,
PrecisionGoal -> 5 CPrec, WorkingPrecision -> 5 CPrec],
5 CPrec]), 5 CPrec];
test = checkval[res];
If[test != NaN,
If[Im[test]/Re[test] > 1*^-5,
Print["Warning: getphi0i resulted in non-negligible Complex \
value at \[Lambda]=", N[\[Lambda], 6], " and A2=", N[A2, 6]]]];
Abs[res]];(*10 \[Mu]m in MeV*)
(*Run entire cell*)
Clear[ResIncrease, s, r, z, domb];
Clear[\[Rho]V, \[Rho]M]
\[Rho]V :=
2.282744303895228569942873684140791998747198853011172394898362081328\
00265957484953105449676513671875`1000.*^-20;
rhoV = \[Rho]V;
\[Rho]M = SetPrecision[2514.0`*kgm32MeV4, 5 CPrec];
rhoM = \[Rho]M;
ResIncrease = 10^0;
s = 1/ResIncrease;
r[z_] := m Tan[s*z];
z[r_] := ArcTan[r/m]/s;
m = m2invMeV;
dd = SetPrecision[ArcTan[10.0`*^-6]/s, 5 CPrec]; (* plate separation *)
dp = SetPrecision[ArcTan[100.0`*^-6]/s,
5 CPrec]; (*thickness of upper plate*)
ddr = 10.0`*^-6*m2invMeV;
dpr = 100.0`*^-6*m2invMeV
initialval1[r_,
A2_?NumberQ, \[Lambda]_?NumberQ, \[Beta]_?
NumberQ] := \[Phi]2[\[Lambda],
A2, \[Beta], \[Phi]0[\[Lambda], A2, \[Beta],
ddr/2, \[Rho]M, \[Rho]V], ddr/2, r, \[Rho]M, \[Rho]V]
Seed[z_, A2_?NumberQ, \[Lambda]_?NumberQ, \[Beta]_?NumberQ] :=
initialval1[r[z], A2, \[Lambda], \[Beta]];
Clear[rho];
rho[z_] :=
SetPrecision[
Piecewise[{{\[Rho]V, (-0.5*dd < z < 0.5*dd)}, {\[Rho]M,
z <= -0.5*dd || z >= -0.5*dd}}], 5 CPrec];
\[Rho][z_] := rho[z]
Clear[C1, C2, C3, C4];
C1[z_?NumberQ] := SetPrecision[mpl/(m^2*s^2) Cos[s*z]^4, 5 CPrec];
C2[z_?NumberQ] :=
SetPrecision[-2 mpl/(m^2*s) Cos[s*z]^4 Tan[s*z], 5 CPrec];
C3[z_?NumberQ, A2_?NumberQ] :=
SetPrecision[-A2/mpl \[Rho][z], 5 CPrec];
C4[f_, \[Lambda]_?NumberQ, A2_?NumberQ, \[Beta]_?NumberQ] :=
If[- ((\[Lambda] f)/mpl) + Log[\[Lambda]] + Log[10]*\[Beta] < 0,
SetPrecision[E^
SetPrecision[-Sqrt[(- ((\[Lambda] f)/mpl) + Log[\[Lambda]] +
Log[10]*\[Beta])^2], 5 CPrec], 5 CPrec],
SetPrecision[E^
SetPrecision[
Sqrt[(- ((\[Lambda] f)/mpl) + Log[\[Lambda]] +
Log[10]*\[Beta])^2], 5 CPrec], 5 CPrec]]
Residual[z_, A2_, \[Lambda]_, \[Beta]_, u_] :=
Abs[C1[z]*u''[z] + C2[z]*u'[z] + C3[z, A2]*u[z] +
C4[u[z], \[Lambda], A2, \[Beta]]]
RelResidual[z_, A2_, \[Lambda]_, \[Beta]_, u_] :=
Residual[z, A2, \[Lambda], \[Beta], u]/
Max[Abs[C1[z]*u''[z]], Abs[C2[z]*u'[z]], Abs[C3[z, A2]*u[z]],
Abs[C4[u[z], \[Lambda], A2, \[Beta]]]]
coordlist[dmin_, ddelta_, dfactor_, nmin_, nmax_] :=
Block[{i, dummy},
dummy = Table[N[dmin, 5 CPrec], nmax - nmin + 1];
If[nmin == 0,
For[i = 1, i <= nmax, i++,
dummy[[i + 1]] = dummy[[i]] + N[ddelta*dfactor^i, 5 CPrec];
];,
dummy[[1]] += Sum[N[ddelta*dfactor^i, 5 CPrec], {i, 1, nmin}];
For[i = 2, i <= nmax - nmin + 1, i++,
dummy[[i]] =
dummy[[i - 1]] + N[ddelta*dfactor^(nmin + i - 1), 5 CPrec];
];];
N[dummy, 5 CPrec]
];
Clear[bvp, x, \[Lambda], A2, \[Beta]];
bvp = {SetPrecision[
C1[z]*D[D[u[z], z], z] + C2[z]*D[u[z], z] + C3[z, A2]*u[z] +
C4[u[z], \[Lambda], A2, \[Beta]], 5 CPrec] ==
SetPrecision[0, 5 CPrec],
SetPrecision[u[SetPrecision[-\[Pi]/2, 5 CPrec]], 5 CPrec] ==
N[SetPrecision[\[Phi]M[A2, \[Lambda], \[Beta], \[Rho]M],
5 CPrec], 5 CPrec],
SetPrecision[u[SetPrecision[\[Pi]/(2*s), 5 CPrec]], 5 CPrec] ==
SetPrecision[
N[\[Phi]M[A2, \[Lambda], \[Beta], \[Rho]M], 5 CPrec],
5 CPrec]} /. m -> m2invMeV;
dist1 = SetPrecision[\[Pi]/(2.0`*s),
5 CPrec];(* from -\[Infinity] to the border of the lower plate*)
(* Was ist dist 2? dp/2 nicht legitim?*)
dist2 = SetPrecision[dp/2., 5 CPrec];(* half thickness of the plate*)
(*was sollen die kommenden Größen sein?*)
nMaxRes = 110*^5;(*max. resolution in the gap 110*^6*)
npts1 = 80;(*points in free vacuum and free material, 400*)
npts2 = 80;(*points in the half-gap, 500*)
npts3 = 80;(*points in the upper plate, 500*)
nFinePoints = 30;(*120*)
minMeshSize = SetPrecision[(s*dd/nMaxRes), 5 CPrec];
f1 =(*von -\[Pi]/2 bis -dd/2*)
SetPrecision[(\[Delta] /.
Flatten[FindRoot[
SetPrecision[
Sum[N[minMeshSize, 5 CPrec]*\[Delta]^N[n, 5 CPrec], {n, 1,
npts1 - 1}], 5 CPrec] ==
SetPrecision[
dist1 - dd/2 - N[nFinePoints*minMeshSize, 5 CPrec],
5 CPrec], {\[Delta], 1.01`}, MaxIterations -> 3000,
PrecisionGoal -> 5 CPrec + 2, AccuracyGoal -> 5 CPrec + 2,
WorkingPrecision -> 5 CPrec]][[1]]), 5 CPrec];
(*f2= (*half thickness of the upper \
plate*)SetPrecision[(\[Delta]/.Flatten[FindRoot[SetPrecision[Sum[N[\
minMeshSize,5 CPrec]*\[Delta]^N[n,5 CPrec],{n,1,npts3-1}],5 CPrec]\
\[Equal]SetPrecision[dist2-N[nFinePoints*minMeshSize,5 CPrec],5 \
CPrec],{\[Delta],1.01`},MaxIterations\[Rule]3000,PrecisionGoal\[Rule]\
5 CPrec+2,AccuracyGoal\[Rule]5 CPrec+2,WorkingPrecision\[Rule]5 \
CPrec]][[1]]),5CPrec];*)
(*f3=(*points above upper plate divided by \
2*)SetPrecision[(\[Delta]/.Flatten[FindRoot[SetPrecision[Sum[N[\
minMeshSize,5 CPrec]*\[Delta]^N[n,5 CPrec],{n,1,npts1-1}],5 CPrec]\
\[Equal]SetPrecision[0.5*(dist1-dp-dd/2)-N[nFinePoints*minMeshSize,5 \
CPrec],5 \
CPrec],{\[Delta],1.01`},MaxIterations\[Rule]3000,PrecisionGoal\[Rule]\
5 CPrec+2,AccuracyGoal\[Rule]5 CPrec+2,WorkingPrecision\[Rule]5 \
CPrec]][[1]]),5CPrec];*)
f4 =(*half gab between plates*)
SetPrecision[(\[Delta] /.
Flatten[FindRoot[
SetPrecision[
Sum[N[minMeshSize, 5 CPrec]*\[Delta]^N[n, 5 CPrec], {n, 1,
npts2 - 1}], 5 CPrec] ==
SetPrecision[dd/2 - N[nFinePoints*minMeshSize, 5 CPrec],
5 CPrec], {\[Delta], 1.001`}, MaxIterations -> 3000,
PrecisionGoal -> 5 CPrec + 2, AccuracyGoal -> 5 CPrec + 2,
WorkingPrecision -> 5 CPrec]][[1]]), 5 CPrec];
points = SetPrecision[Sort[N[Join[
(* lower plate *)
coordlist[
N[-nFinePoints*minMeshSize - dd/2, 5 CPrec], -minMeshSize, f1,
1, npts1 - 1],
Table[
N[-i*minMeshSize - dd/2, 5 CPrec], {i, 0, nFinePoints - 1}],
Table[N[i*minMeshSize - dd/2, 5 CPrec], {i, 1, nFinePoints}],
Table[-\[Pi]/2 + N[i*minMeshSize, 5 CPrec], {i, 1,
nFinePoints}] ,
(*in between plates*)
coordlist[N[nFinePoints*minMeshSize - dd/2, 5 CPrec],
minMeshSize, f4, 1, npts2 - 2],
coordlist[
dd/2 - N[nFinePoints*minMeshSize, 5 CPrec], -minMeshSize, f4,
1, npts2 - 1],
Table[
dd/2 - N[i*minMeshSize, 5 CPrec], {i, 0, nFinePoints - 1}](*
inter-spacing *),
Table[dd/2 + N[i*minMeshSize, 5 CPrec], {i, 1, nFinePoints - 1}],
(* upper plate *)
coordlist[N[nFinePoints*minMeshSize + dd/2, 5 CPrec],
minMeshSize, f1, 1, npts1 - 1],
Table[\[Pi]/2 - N[i*minMeshSize, 5 CPrec], {i, 1, nFinePoints}]
]]](*free space from right*), 5 CPrec];
points[[1]] = -\[Pi]/2;
points[[-1]] = \[Pi]/(2*s);
domb = SetPrecision[
ToElementMesh["Coordinates" -> Transpose[{points}],
"MeshElements" -> {LineElement[
Table[{i, i + 1}, {i, 1, Length[points] - 1}]]},
"BoundaryElements" -> {(*external boundaries:*)
PointElement[{{1}, {Length[points]}}]}], 5 CPrec];
Length[points]
After the above code has been activated one can simply solve the ODE by deciding on the parameters (following code). This gives the correct result for [Lambda]=10^(27), not for [Lambda]=10^(31)
Clear[\[Lambda]t, A2, \[Lambda], A2, Y1, Y2, \[Beta]]
\[Lambda] = 10^31;
A2 = 10^40;
\[Beta] = 1;
Y1 = NDSolveValue[Evaluate[bvp], u, z \[Element] domb,
InitialSeeding -> {u[z] == Seed[z, A2, \[Lambda], \[Beta]]}]
Plot[{Y1[z], Seed[z, A2, \[Lambda], \[Beta]]}, {z, -dd, dd}]
Any help or insight would be appreciated, thank you!
