In[2993]:= 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 = 7; r1 = .2; rgamma = 1; rbeta = .6; k = 1; xbeta = .4; \
beta = .2; sbar = .2; zeta = 5;
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] 

During evaluation of In[2993]:= b = 0.965756

During evaluation of In[2993]:= w = 0.842599

During evaluation of In[2993]:= theta = 0.979777

Out[3007]= Null^3

During evaluation of In[2993]:= sstar = 0.2

During evaluation of In[2993]:= pstar = 0.590558

Out[3008]= Null^2

During evaluation of In[2993]:= welfare=0.965116

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