In[3028]:= Clear["Global`*"]
theta = (h - p*h + p*b)/h ;
 
 gamma = 1 - beta;

w = (gamma * rgamma + theta* beta* (rbeta - xbeta))/ (gamma + 
     theta *beta) ;

 

welfare = 
  gamma*r1 + 
   gamma*w + (b/h)*beta *r1 + (b/h)*beta *
    w + (h^2*beta)/(2*h) - (b^2*beta)/(2*h) +
   (1 - b/h)*
    beta * ((1 - p)*r1 + p*(1 - s)*r1 - p*s*d + (1 - p)*w + 
      p * rbeta) - (1 - b/h)*beta * h - (1 - b/h)*beta *p*s*k   - 
   p^2/zeta ;

 

b = Min [1, 
   1/(2 beta) (-h + beta h p + beta d p s + beta p r1 s - 
      beta p xbeta + \[Sqrt]((h - beta h p - beta d p s - 
           beta p r1 s + beta p xbeta)^2 - 
         4 beta (h p rbeta - beta h p rbeta - h p rgamma + 
            beta h p rgamma - d h p s + beta d h p^2 s - h p r1 s + 
            beta h p^2 r1 s + beta h p xbeta - beta h p^2 xbeta)))];
h = 1; d = 1; r1 = .2; rgamma = 1; rbeta = .6; k = 20; xbeta = .4; \
beta = .5; sbar = .2; zeta = 5;
temp = NMaximize[{welfare, s > 0, s < sbar, p > 0, p < 1}, {s, p}]
temp = NMaximize[{welfare, s > 0, s < sbar, p > 0, p < 1}, {s, p}] ;  
 (*p=.3 ; s=.3;*)
(*Show[ContourPlot[welfare,{s,0, .5}, {p,0, 1}  ]]*)
temp1 = Part[temp, 2];
temp2 = Part[temp1, 1];     
   sstar = s /. %;
temp3 = Part[temp1, 2];
pstar = p /. %;
 s = sstar; p = pstar;
Print["b = ", b]   Print["w = ", w] Print ["theta = ", theta]
Print["sstar = ", sstar] Print["pstar = ", pstar] 
Print["welfare=", welfare] 

Out[3035]= {0.630277, {s -> 0., p -> 0.704317}}

During evaluation of In[3028]:= b = 0.117099

During evaluation of In[3028]:= w = 0.780485

During evaluation of In[3028]:= theta = 0.378158

Out[3043]= Null^3

During evaluation of In[3028]:= sstar = 0.

During evaluation of In[3028]:= pstar = 0.704317

Out[3044]= Null^2

During evaluation of In[3028]:= welfare=0.630277

Solve[ b == p *(w - rbeta) + p*s*(r1 + d), b]