Hello

Here is the case I don't know how to solve: Over a range of "R", I want to solve etha=Exp[function of etha and R]. That should give an "etha" for each "R" Then calculate some other parameters that are function of both "etha" and "R" I really appreciate any help The code is

Clear["Global`*"]
k = 1.3806488*10^-23;      
Av = 6.02*10^23;                   
T = 25 + 273.15;                       
Psat1 = 3.168*10^3;          
PL = 101.325*10^3;             
vLinf = (1/997.0479)*18.015/1000;   
gama = 0.07199;                      
theta = Degree 10 ;                
qc = 1.22*10^8;
W = 50*10^-6;                        
alfa = Degree 45;
L = 0.45;                                    
Vcontainer = L^3 ;             
N1 = 3.049*10^27;             
N2 = 3.67*10^22;               
KH = 9.077*10^4;        
bound1 = W/Cos [theta - alfa];
Print["bound1", "=", N[bound1]]
(***********************************************************)
Vv[R_] := (Pi*R^3/3)*(2 - 
    3*Sin[theta - alfa] + (Sin[
       theta - alfa])^3 + (Cos[theta - alfa])^3/Tan[alfa])
Asv[R_] := Pi*R^2*(Cos[theta - alfa])^2/Sin[alfa];
Alv[R_] := 2*Pi*R^2*(1 - Sin[theta - alfa]);
H[R_] := R*(1 - Sin[theta - alfa] + Cos[theta - alfa]/Tan[alfa]);

etha[R_] := 
  Exp[(vLinf/(6.02*10^23))*(PL - Psat1)/(k*T) - 
    N2*k*T/(N1*k*T - qc*Vv[R]*(etha[R]*Psat1 - KH))];

N1L[R_] := N1 - qc*(etha[R]*Psat1*Vv[R]/(k*T));     
N2sL[R_] := (N1L[R]*PL/KH);    
N2L[R_] := N2/(1 + qc*PL*Vv[R]/(N2sL[R]*k*T));    

Rc[R_] := 2*gama/(etha*Psat1 + PL*(N2L[R]/N2sL[R]) - PL);    
B[R_] := qc*(-2*gama * Vv[R]/Rc[R] + 
      gama*(Asv[R]* Cos [theta] + Alv[R])) + 
   N1*k*T*(N2/N1 - N2L[R]/N1L[R]) + 
   N2*k*T*Log[N1*k*T/(N1*k*T - qc*etha*Psat1*Vv[R] + qc*Vv[R]*KH)];


(* I want to plot B vs. H same as Fig 2*)
t2 = Table[{ H[R], B[R]/qc}, {R, 10^-6, bound1, 10^-7}];
ListPlot[t2, Joined -> True, PlotStyle -> {Thick}, 
 AxesLabel -> {Style["H", "Graphics", FontSize -> 13], 
   Style["B-B0 /qc (J)", "Graphics", FontSize -> 13]}, 
 AxesOrigin -> {0, 0}, TicksStyle -> Directive[12]]