I'm more accustomed to working in a slightly different format, but this should be equivalent to what you're doing. Define functions to build the pdes and g:

pdef[\[Beta]_, \[Gamma]_, \[Sigma]_, \[Mu]_, \[Rho]_, \[Zeta]_] := 
  2 \[Beta] f[r] + 
    2 r (-\[Mu] + \[Rho]) f'[
      r] + (2 (-1)^(\[Gamma]/(-1 + \[Gamma])) r^((2 \[Gamma])/(-1 + \
\[Gamma])) (-1 + \[Gamma]) \[Zeta]^(-(\[Gamma]/(-1 + \[Gamma]))) \
\[Rho]^(\[Gamma]/(-1 + \[Gamma])) f'[
         r]^(\[Gamma]/(-1 + \[Gamma])))/\[Gamma] - 
    r^2 \[Sigma]^2 f''[r] == 0;

g[\[Beta]_, \[Gamma]_, \[Mu]_, \[Eta]_, \[Rho]_, \[Sigma]_, \
\[Lambda]_][
   r_] := (2^(1 - \[Gamma]) (1/
        r)^\[Gamma] \[Eta]^\[Gamma] \[Rho]^-\[Gamma] ((-2 \[Beta] + \
\[Gamma] (-2 \[Mu] + 
             2 \[Rho] + (1 + \[Gamma]) \[Sigma]^2))/(-1 + \
\[Gamma]))^(-1 + \[Gamma]))/\[Gamma];

The solver with most of the constants plugged in. (You don't seem to use J or k for anything as far as I can tell.)

solver[rSol_?NumericQ] := 
  NDSolve[{pdef[0.012, 1/2, 0.2, 0.03, 0.02, 1/2], 
    f[rSol] == g[0.012, 1/2, 0.03, 0.8, 0.02, 0.2, 1.03][rSol], 
    Derivative[1][f][rSol] == 
     Derivative[1][g[0.012, 1/2, 0.03, 0.8, 0.02, 0.2, 1.03]][rSol]}, 
   f, {r, rSol/10, rSol}];

The main difference is that I removed the Dirichlet code and just used a simple f[rSol] = g[..][rSol] and the corresponding derivatives directly, using the conditions you wrote in the text.

FindRoot complained about the solution being non-numeric, so I created an auxillary function to keep it happy:

rsf[rs_?NumericQ] := 
  Derivative[1][f][1/10^5] /. Flatten@solver[1/10^5];

I changed the constant to a value different from yours, to verify that it did indeed search for the root:

FindRoot[
 rsf[rSol] - 
  Derivative[1][g[0.012, 1/2, 0.03, 0.8, 0.02, 0.2, 1.03]][
   rSol], {rSol, .0002(*1/100000*)}]

Out[39]= {rSol -> 0.00001}

I reran once last time to generate plots and verify other conditions. The output was of the form {f->InterpolatingFunction[..]}, which I won't copy here since it's long and ugly.

nds = Flatten@solver[.00001]

Here are f,g, and f-g:

Since the scales are quite different, here's f-g alone:

Here's f'[r]

[

And here are {{f[rSol],g[rSol]}{f'[rSol],g'[rSol]}}. f,g and f',g' are equal at rSol, as required.

{f[r] /. nds, 
   g[0.012, 1/2, 0.03, 0.8, 0.02, 0.2, 1.03][r]}, {f'[r] /. nds, 
   D[g[0.012, 1/2, 0.03, 0.8, 0.02, 0.2, 1.03][r], r]}} /. 
 r -> 9.999999999999999`*^-6

Out[53]= {{63245.6, 63245.6}, {-3.16228*10^9, -3.16228*10^9}}