Clear[fdd, pdetoode, tooderule, sollst]
fdd[{}, grid_, value_, order_] := value;
fdd[a__] := NDSolve`FiniteDifferenceDerivative@a;
pdetoode[funcvalue_List, rest__] :=
pdetoode[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[[1]],
rest];
pdetoode[{func__}[var__], rest__] :=
pdetoode[Alternatives[func][var], rest];
pdetoode[rest__, grid_?VectorQ, o_Integer] :=
pdetoode[rest, {grid}, o];
pdetoode[func_[var__], time_, {grid : {__} ..}, o_Integer] :=
With[{pos = Position[{var}, time][[1, 1]]},
With[{bound = #[[{1, -1}]] & /@ {grid},
pat = Repeated[_, {pos - 1}],
spacevar = Alternatives @@ Delete[{var}, pos]},
With[{coordtoindex =
Function[coord,
MapThread[
Piecewise[{{1, # === #2[[1]]}, {-1, # === #2[[-1]]}},
All] &, {coord, bound}]]},
tooderule@
Flatten@{((u : func) |
Derivative[dx1 : pat, dt_, dx2___][(u : func)])[x1 : pat,
t_, x2___] :> (Sow@coordtoindex@{x1, x2};
fdd[{dx1, dx2}, {grid},
Outer[Derivative[dt][u@##]@t &, grid],
"DifferenceOrder" -> o]),
inde : spacevar :>
With[{i = Position[spacevar, inde][[1, 1]]},
Outer[Slot@i &, grid]]}]]];
tooderule[rule_][pde_List] := tooderule[rule] /@ pde;
tooderule[rule_]@Equal[a_, b_] :=
Equal[tooderule[rule][a - b], 0] //.
eqn : HoldPattern@Equal[_, _] :> Thread@eqn;
tooderule[rule_][expr_] := #[[Sequence @@ #2[[1, 1]]]] & @@
Reap[expr /. rule]
Clear@pdetoae;
pdetoae[funcvalue_List, rest__] :=
pdetoae[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[[1]], rest];
pdetoae[{func__}[var__], rest__] :=
pdetoae[Alternatives[func][var], rest];
pdetoae[func_[var__], rest__] :=
Module[{t},
Function[pde, #[
pde /. {Derivative[d__][u : func][inde__] :>
Derivative[d, 0][u][inde, t], (u : func)[inde__] :>
u[inde, t]}] /. (u : func)[i__][t] :> u[i]] &@
pdetoode[func[var, t], t, rest]]
Clear@k;
Clear@y;
Clear@z;
Clear@\[Phi];
inf = 1000;
eqn1 = y''[x] +
y'[x]/x - ((k - z[x])^2* y[x])/x^2 + ((1 - y[x]^2)* y[x])/2 == 0;
eqn2 = z''[x] - z'[x]/x + (k - z[x])*y[x]^2 == 0;
eqns = Join[{eqn1, eqn2}];
bc = {y[0] == 0, y[inf] == 1, z[0] == 0, z[inf] == k};
set = Join[eqns, bc];
points = 1000;
domain = {0, inf};
grid = Array[# &, points, domain];
difforder = 4;
(*Definition of pdetoode is above.*)
ptoafunc = pdetoae[{y[x], z[x]}, grid, difforder];
delete = #[[2 ;; -2]] &;
ae = delete /@ ptoafunc[Map[x^2 # &, eqns, {2}] // Simplify];
WorkingPrecision -> 30;
sollst = Table[
ListInterpolation[#, grid] & /@ Partition[#, points] &@
With[{initial = 1},
FindRoot[{ae, bc},
Flatten[{{y@#, initial} & /@ grid, {z@#, initial} & /@ grid},
1]]][[All, -1]], {k, 3}]
Attachments: