The method I am using here is called Successive Over Relaxation (SOR). There are plenty of references that you can find on the web on this. My implementation is a very simple version. There are much more sophisticated methods in SOR that could be used.

Thank you for response. A further extension to this that I have been attempting is to produce plots of the magnetic field, i.e. plotting the following function based on the numerical solutions: bfield = f'[x]/x;

My attempt so far has yielded interesting results that seem to break down at zero.

bfield = f'[x]/x;
bfieldplot = 
 Plot[Table[bfield, {f, zsol}] // Evaluate, {x, 0, 10}, 
  PlotRange -> {0, 0.6}, 
  PlotStyle -> 
   Table[Blend[{Blue, Cyan}, x], {x, 0, 1, 1/(Length[zsol] - 1)}], 
  GridLines -> Automatic, AxesLabel -> {r, B}]

Is there any way to rectify this?

I have also attached my current Mathematica notebook for further reference.

Attachments:

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:

Unfortunately, the numerical solutions associated with your method are actually worse than my original finite difference method which I have shown above, rather than better.