0
$\begingroup$

I can't use my Mathematica 9 code in Mathematica 11 — it works in Mathematica 9, but when I run it in Mathematica 11 it shows a lot of errors It gives me this message

ReplaceAll::reps: {NDSolve[{(2 (ψ^′)[r])/r+(ψ^′′)[r]==0.00135103 Sinh[ψ[r]],ψ[1]==-5.,(ψ^′)[1.4938]==0.001},ψ,{r,1,1.4938},Method->{EquationSimplification->Residual}]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing..

Please help me fix my code or tell me why NDSolve did not work in Mathematica 11?

Code

Clear["Global`*"];
(*initialize variables*)
Subscript[λ, B] = 0.714107`20;  
 a = 25;  
Z = 1000;  
Subscript[c, 0] = .2; 
η = 0.3;  
γ = Subscript[λ, B]/a;  
R = η^(-1/3);  
Subscript[n, m] = 3 η/(4 π);
Avogadro = 6.0221413`20 10^-7; 
Subscript[n, 0] = Subscript[c, 0] Avogadro a^3;  
κ2 = 8 π γ Subscript[n, 0] ; 
κ = Sqrt[κ2] ; 
ϵ = $MachineEpsilon ;

ψin[ψtest_?NumericQ] := 
 NDSolve[{ψ''[r] + 2 ψ'[r]/r == κ2 Sinh[ψ[r]] + 
     3 Z γ, ψ[1] == ψtest, ψ'[ϵ] == 
    0}, ψ, {r, ϵ, 1}, 
  Method -> {"EquationSimplification" -> "Residual"}]

ψout[ψtest_?NumericQ] := 
 NDSolve[{ψ''[r] + 
     2 ψ'[r]/r == κ2 Sinh[ψ[r]], ψ[
     1] == ψtest, ψ'[R] == 0}, ψ, {r, 1, R}, 
  Method -> {"EquationSimplification" -> "Residual"}]

ψinTry[ψtest_?
   NumericQ] := (Evaluate[ψ'[1] /. ψin[ψtest]])[[1]]

ψoutTry[ψtest_?
   NumericQ] := (Evaluate[ψ'[1] /. ψout[ψtest]])[[1]]

ψtest0 = -5;  
ψa = ψtest /. 
  FindRoot[ψinTry[ψtest] == ψoutTry[ψtest], {\
ψtest, ψtest0}];
Clear[ψtest];
ψtest0 = ψa
ψR = (Evaluate[ψ[R] /. ψout[ψa]])[[1]];
Subscript[c, s] = 
 3 Subscript[c, 
  0] η NIntegrate[(r^2) If[
     r < 1, (Evaluate[Exp[ψ[r] /. ψin[ψa]]])[[
      1]], (Evaluate[Exp[ψ[r] /. ψout[ψa]]])[[1]]], {r,
     0, R}] 

Print[NumberForm[η, 4], "  ", NumberForm[Subscript[c, s], 10], 
  "  ", NumberForm[ψR, 10]];

Clear[ψtest];
Clear[r];

yin = ψin[ψa];
yout = ψout[ψa];

ψlist = {};
Elist = {};

ψr[r_] := If[r < 1, ψ[r] /. yin[[1]], ψ[r] /. yout[[1]]];

Plot[ψr[r], {r, 0, R}, PlotRange -> All, 
 AxesLabel -> {"r/a", "ψ(r)"}, 
 BaseStyle -> {FontWeight -> "Bold", FontFamily -> "Calibri", 
   FontSize -> 16}, 
 EvaluationMonitor :> AppendTo[ψlist, {r, ψr[r]}]]

Export["psi.dat", Sort[ψlist], "Table"]

Er[r_] := If[r < 1, ψ'[r] /. yin[[1]], ψ'[r] /. yout[[1]]];

Plot[Er[r], {r, 0, R}, PlotRange -> All, AxesLabel -> {"r/a", "E(r)"},
  BaseStyle -> {FontWeight -> "Bold", FontFamily -> "Calibri", 
   FontSize -> 16}, EvaluationMonitor :> AppendTo[Elist, {r, Er[r]}]]

Export["E.dat", Sort[Elist], "Table"]
$\endgroup$
  • $\begingroup$ \[Psi]test is not defined in the line FindRoot[\[Psi]inTry[\[Psi]test]==\[Psi]outTry[\[Psi]test],{\[Psi]test,\[Psi]test0}] so I have no idea how this could have worked in any Mathematica version. All what you had to do is run each command in separate cell, until you see where the problem is. $\endgroup$ – Nasser Dec 3 '17 at 0:16
  • $\begingroup$ You defined \[Psi]inTry[\[Psi]test_?NumericQ] as taking numeric. But you call it using \[Psi]inTry[\[Psi]test] where \[Psi]test is not defined. Same for \[Psi]outTry[\[Psi]test_?NumericQ]. Your code seems to be missing definitions. $\endgroup$ – Nasser Dec 3 '17 at 0:34
  • $\begingroup$ Thank you for your help $\endgroup$ – user53936 Dec 3 '17 at 21:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.