This problem can be solve with FEM and FDM as well. First we show FDM code since pdetoode is implementation of FDM. We copy this function as well from here (thanks to xzczd)

Clear[fdd, pdetoode, tooderule, pdetoae, diffbc, rebuild]
fdd[{}, grid_, value_, order_, periodic_] := 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[front__, grid_?VectorQ, o_Integer, periodic_ : False] := 
  pdetoode[front, {grid}, o, periodic];

pdetoode[func_[var__], time_, {grid : {__} ..}, o_Integer, 
   periodic : True | False | {(True | False) ..} : False] := 
  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, PossibleZeroQ[# - #2[[1]]]}, {-1, 
             PossibleZeroQ[# - #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, 
           PeriodicInterpolation -> periodic]), 
        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]

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]]

diffbc[rst__][a : _List | _Equal] := diffbc[rst] /@ a
diffbc[dvar : {t_, order_} | (t_) .., sf_ : 0][a_] /; sf =!= t := 
 sf a + D[a, dvar]

rebuild[funcarray_, grid_?VectorQ, timeposition_ : 1] := 
 rebuild[funcarray, {grid}, timeposition]

rebuild[funcarray_, grid_, timeposition_?Negative] := 
 rebuild[funcarray, grid, Range[Length@grid + 1][[timeposition]]]

rebuild[funcarray_, grid_, timeposition_ : 1] /; 
  Dimensions@funcarray === Length /@ grid := 
 With[{depth = Length@grid}, 
  ListInterpolation[
     Transpose[
      Map[Developer`ToPackedArray@#["ValuesOnGrid"] &, #, {depth}], 
      Append[Delete[Range[depth + 1], timeposition], timeposition]], 
     Insert[grid, Flatten[#][[1]]["Coordinates"][[1]], 
      timeposition]] &@funcarray]

(*Definitions*)(*Constants*)\[Mu] = 
  2.75*10^-6 (*mN mm^-2 s*);(*Viscosity of water*)
\[Gamma] = 0.072(*mN mm^-1*);(*Air-water surface tension*)
\[Beta] = 1(*mm s^-1*);(*Slip constant*)

(*Grid specification*)
difforder = 4;
lb = 1/100; rb = 1; points = 20;
grid = Array[# &, points, {lb, rb}];

unitStepExpand = Simplify`PWToUnitStep@PiecewiseExpand@# &;

With[{h = h[r/R[t], t], q = q[r/R[t], t], p = p[r/R[t], t], 
  R = R[t]},(*Equations neglecting evaporation*)
 EqnP = Simplify[
    p == \[Gamma] (-D[h, r, r] - 
        unitStepExpand@If[x < 1/20, D[h, r, r], D[h, r]/r])] /. 
   r -> x R;
 EqnQ = Simplify[q == (r h^3)/(3 \[Mu]) (-D[p, r])] /. r -> x R;
 EqnC = Simplify[D[h, t] == -(1/r) D[q, r]] /. r -> x R;
 EqnR = D[R, t] == \[Beta] (-D[h, r]^3 - 1) /. r -> R;]

(*Boundary conditions*)
dr = R[t] lb;

EqnBC1 = {D[r h[lb, t] dr, t] + r q[lb, t] == 0} /. r -> lb R[t];
EqnBC2 = {Derivative[1, 0][h][lb, t] == 0};
EqnBC3 = {h[rb, t] == 0};

(*Initial conditions*)
EqnIC = {h[x, 0] == (1 - x^2)};

(*Equation Discretisation*)
removeredundant = #[[3 ;; -2]] &;

tfunc = pdetoode[{p, q, h}[x, t], t, grid, difforder];
odemid = Map[tfunc, {EqnP, EqnQ}, {2}];
odeC = Block[{p, q}, Set @@@ odemid; tfunc@EqnC] // removeredundant;

odeBC1 = Block[{p, q}, Set @@@ odemid; tfunc@EqnBC1];
odeBC2 = With[{sf = 100}, diffbc[t, sf]@EqnBC2 // tfunc];
odeBC3 = With[{sf = 100}, diffbc[t, sf]@EqnBC3 // tfunc];

odeR = Block[{p, q}, Set @@@ odemid; tfunc@EqnR];

odeIC = Append[EqnIC // tfunc, R[0] == 10];

FDM solution of the Laplace equation in general form

delx = Differences[grid][[1]]; 
XYgrid[dom_List, pts_List] := 
 MapThread[
  N@Range[Sequence @@ #1, Abs[Subtract @@ #1]/#2] &, {dom, pts - 1}];
BoundaryIndex[xgridlen_, ygridlen_] := 
  Module[{tmp, left, right, bot, top}, 
   tmp = Table[(n - 1) ygridlen + Range[1, ygridlen], {n, 1, 
      xgridlen}]; {left, right} = tmp[[{1, -1}]]; {bot, top} = 
    Transpose[{First[#], Last[#]} & /@ tmp]; {top, right[[2 ;; -2]], 
    bot, left[[2 ;; -2]]}];
FDMat[deriv_, xygrid_, difforder_] := 
 Map[NDSolve`FiniteDifferenceDerivative[#, xygrid, 
     "DifferenceOrder" -> difforder]["DifferentiationMatrix"] &, deriv]
{domain, pts, 
   difforder} = {{{lb, 5 points delx + lb }, {0, H}}, {5 points + 1, 
    20}, 4};
xygrid = XYgrid[domain, pts]; {nx, ny} = 
 Map[Length, xygrid]; {top, right, bot, left} = 
 BoundaryIndex[nx, ny]; {dx, dy, dx2, dy2} = 
 FDMat[{{1, 0}, {0, 1}, {2, 0}, {0, 2}}, xygrid, 
  difforder]; boundaries = Join[top, right, bot, left]; sgrid = 
 Flatten[Outer[List, Sequence @@ xygrid], 1]; bot1 = 
 Take[bot, points]; bot2 = Complement[bot, bot1];
u = Table[uu[i], {i, nx ny}]; X = 
 sgrid[[All, 1]]; eqs = (dx2 + dy2) . u + (dx . u)/X; 
eqs[[left]] = (dx . u)[[left]]; eqs[[right]] = u[[right]]; 
eqs[[top]] = u[[top]];
varh = h[#][t] & /@ grid; eqs[[bot1]] = u[[bot1]] - varh; 
eqs[[bot2]] = (dy . u)[[bot2]];

{vec, mat} = CoefficientArrays[eqs, u];

inv = Inverse[mat // N];

U = -inv . vec;

ju = (dy . U)[[bot1]]; ju1 = (ju // removeredundant)/R[t];

odeC1 = Thread[odeC[[All, 1]] == ju1];

Please note, that ju is the evaporation rate expressed as a function of h[t] on the grid. We use odeC1 to compute solution as follows

Monitor[sollst = 
   NDSolveValue[{odeC1, odeR, odeBC1, odeBC2, odeBC3, 
     odeIC}, {h /@ grid, R}, {t, 0, 6}, 
    EvaluationMonitor :> (time = t), 
    Method -> {"EquationSimplification" -> "Solve"}], time];


{hsol} = rebuild[#, grid, 2] & /@ {sollst[[1]]};
Rsol = sollst[[2]];

Animation

frames = 
  Table[Plot[hsol[Abs[r]/Rsol[t], t], {r, 1/10, Rsol[t]}, 
    PlotRange -> {{0, 10}, {0, 1}}], {t, 0, 6, .05}];