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 ([Psi]^[Prime])[r])/r+([Psi]^[Prime][Prime])[r]==0.00135103 Sinh[[Psi][r]],[Psi][1]==-5.,([Psi]^[Prime])[1.4938]==0.001},[Psi],{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.-. If you can help me fixing my code. OR let me know why ND-solve did not work in Mathematica 11 ?

## code

`Clear["Global`*"]; (*initialize variables*) Subscript[\[Lambda], B] = 0.714107`20; a = 25; Z = 1000; Subscript[c, 0] = .2; \[Eta] = 0.3; \[Gamma] = Subscript[\[Lambda], B]/a; R = \[Eta]^(-1/3); Subscript[n, m] = 3 \[Eta]/(4 \[Pi]); Avogadro = 6.0221413`20 10^-7; Subscript[n, 0] = Subscript[c, 0] Avogadro a^3; \[Kappa]2 = 8 \[Pi] \[Gamma] Subscript[n, 0] ; \[Kappa] = Sqrt[\[Kappa]2] ; \[Epsilon] = $ MachineEpsilon ; \[Psi]in[\[Psi]test_?NumericQ] := NDSolve[{\[Psi]''[r] + 2 \[Psi]'[r]/r == \[Kappa]2 Sinh[\[Psi][r]] + 3 Z \[Gamma], \[Psi][1] == \[Psi]test, \[Psi]'[\[Epsilon]] == 0}, \[Psi], {r, \[Epsilon], 1}, Method -> {"EquationSimplification" -> "Residual"}] \[Psi]out[\[Psi]test_?NumericQ] := NDSolve[{\[Psi]''[r] + 2 \[Psi]'[r]/r == \[Kappa]2 Sinh[\[Psi][r]], \[Psi][ 1] == \[Psi]test, \[Psi]'[R] == 0}, \[Psi], {r, 1, R}, Method -> {"EquationSimplification" -> "Residual"}] \[Psi]inTry[\[Psi]test_? NumericQ] := (Evaluate[\[Psi]'[1] /. \[Psi]in[\[Psi]test]])[[1]] \[Psi]outTry[\[Psi]test_? NumericQ] := (Evaluate[\[Psi]'[1] /. \[Psi]out[\[Psi]test]])[[1]] \[Psi]test0 = -5; \[Psi]a = \[Psi]test /. FindRoot[\[Psi]inTry[\[Psi]test] == \[Psi]outTry[\[Psi]test], {\ \[Psi]test, \[Psi]test0}]; Clear[\[Psi]test]; \[Psi]test0 = \[Psi]a \[Psi]R = (Evaluate[\[Psi][R] /. \[Psi]out[\[Psi]a]])[[1]]; Subscript[c, s] = 3 Subscript[c, 0] \[Eta] NIntegrate[(r^2) If[ r < 1, (Evaluate[Exp[\[Psi][r] /. \[Psi]in[\[Psi]a]]])[[ 1]], (Evaluate[Exp[\[Psi][r] /. \[Psi]out[\[Psi]a]]])[[1]]], {r, 0, R}] Print[NumberForm[\[Eta], 4], " ", NumberForm[Subscript[c, s], 10], " ", NumberForm[\[Psi]R, 10]]; Clear[\[Psi]test]; Clear[r]; yin = \[Psi]in[\[Psi]a]; yout = \[Psi]out[\[Psi]a]; \[Psi]list = {}; Elist = {}; \[Psi]r[r_] := If[r < 1, \[Psi][r] /. yin[[1]], \[Psi][r] /. yout[[1]]]; Plot[\[Psi]r[r], {r, 0, R}, PlotRange -> All, AxesLabel -> {"r/a", "\[Psi](r)"}, BaseStyle -> {FontWeight -> "Bold", FontFamily -> "Calibri", FontSize -> 16}, EvaluationMonitor :> AppendTo[\[Psi]list, {r, \[Psi]r[r]}]] Export["psi.dat", Sort[\[Psi]list], "Table"] Er[r_] := If[r < 1, \[Psi]'[r] /. yin[[1]], \[Psi]'[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"] `