Here's a polynomial interpolation method, which can be be found in Chapter 5 of [Boyd (2014)](http://epubs.siam.org/doi/book/10.1137/1.9781611973525).


    nn = 64;
    z0 = w1 + w2;
    rr = 1.1 w1;
    ff = N[WeierstrassP[z0 + rr #, inv] - L, Precision[#]] &;
    wprec = MachinePrecision;
    tj = 2 Pi*Range[0, nn - 1]/nn;
    wj = N[Exp[I tj], wprec];
    fj = ff /@ wj; (* f[zj] *)
    aa = InverseFourier[fj]/Sqrt[nn];
    
    (* Rough check of accuracy of interpolation *)
    "condition"@# -> Log10@Ratios[#] &@N@MinMax@Abs@fj
    ip = FromDigits[Reverse@aa, (z - z0)/rr];
    Max@Table[(ff[z] - ip)/ff[z] /. z -> z0 + r Exp[I t] // N // Abs,
      {r, (1/8 - 0.01) rr, rr, rr/8}, {t, 0., 2 Pi, 0.2/r}]
    (*
      "condition"[{0.0694093, 0.75655}] -> {1.03742}
      3.98344*10^-11
    *)
    
    z2 = Eigenvalues@companionMatrix[aa];
    roots = z0 + rr*Select[z2, Abs[#] < 0.999 &]
    (*  {0.776416 + 2.35619 I, 2.36518 + 2.35619 I}  *)
    
    Graphics[{
      EdgeForm@Gray, LightBlue, Disk[ReIm@z0, rr], Gray, Point[ReIm@z0],
      First@plt,
      Black, PointSize[Medium], Point@ReIm[roots]},
     Frame -> True]

![Mathematica graphics](https://i.sstatic.net/m59Gn.png)

Auxiliary code:

For the companion matrix, you can use

    companionMatrix[coeffs_] := 
      Join[SparseArray[
        Band[{1, 2}] -> 1, {Length@coeffs - 2, Length@coeffs - 1}],
       {-coeffs[[;; -2]]/coeffs[[-1]]}
       ];

Or

    companionMatrix = NRoots`CompanionMatrix