Huang,

You're right on that point. The plot should come out OK in spite of the error in the integration. My comment was too hasty. René

Hi Rene,

Thank you for your help. It was indeed a mistake.

But what I was trying to plot is a function that is well defined. I suspect that the definition of jointd1 should not affect the Plot3D of jointf1, or it is not? Thanks.

Huang

abarS := (-2 thB + 
     2 thB^2 xi + \[Sqrt]((2 thB - 2 thB^2 xi)^2 - 
        4 (-thB + thB^2 xi) (2 - 2 thB - xi - 2 thB xi + 3 thB^2 xi + 
           2 thB xi^2 - 2 thB^2 xi^2 + 2 thB Log[thB] - 
           4 thB xi Log[thB] - 2 thB^2 xi Log[thB] + 
           2 thB xi^2 Log[thB] + 4 thB^2 xi^2 Log[thB] - 
           2 thB^2 xi^3 Log[thB] + 2 thB Log[1 - xi] - 
           4 thB xi Log[1 - xi] - 2 thB^2 xi Log[1 - xi] + 
           2 thB xi^2 Log[1 - xi] + 4 thB^2 xi^2 Log[1 - xi] - 
           2 thB^2 xi^3 Log[1 - xi] - 2 thB Log[1 - thB xi] + 
           4 thB xi Log[1 - thB xi] + 2 thB^2 xi Log[1 - thB xi] - 
           2 thB xi^2 Log[1 - thB xi] - 
           4 thB^2 xi^2 Log[1 - thB xi] + 
           2 thB^2 xi^3 Log[1 - thB xi])))/(2 (2 - 2 thB - xi - 
       2 thB xi + 3 thB^2 xi + 2 thB xi^2 - 2 thB^2 xi^2 + 
       2 thB Log[thB] - 4 thB xi Log[thB] - 2 thB^2 xi Log[thB] + 
       2 thB xi^2 Log[thB] + 4 thB^2 xi^2 Log[thB] - 
       2 thB^2 xi^3 Log[thB] + 2 thB Log[1 - xi] - 
       4 thB xi Log[1 - xi] - 2 thB^2 xi Log[1 - xi] + 
       2 thB xi^2 Log[1 - xi] + 4 thB^2 xi^2 Log[1 - xi] - 
       2 thB^2 xi^3 Log[1 - xi] - 2 thB Log[1 - thB xi] + 
       4 thB xi Log[1 - thB xi] + 2 thB^2 xi Log[1 - thB xi] - 
       2 thB xi^2 Log[1 - thB xi] - 4 thB^2 xi^2 Log[1 - thB xi] + 
       2 thB^2 xi^3 Log[1 - thB xi]));

xbar := abarS*xi*(1 - thB)/(1 - xi*thB);
c := 1/(abarS*xi*(1 - thB) + xi*thB - 
     abarS*thB*xi*(1 - xi)*Log (thB*(1 - xi)/(1 - xi*thB)));
f1[x_, thB_, xi_] := 
 Piecewise[{{c*(1 - abarS*thB*xi*(1 - xi)/(abarS*xi - x)), 
    0 <= x <= xbar}, {xbar <= x <= 1, c*xi*thB}}, 0]
jointf1[x1_, x2_, thB_, xi_] := 
  Piecewise[{{1 >= x1 > x2, 2*f1[x1, thB, xi]*f1[x2, thB, xi]}}, 0];
jointd1[a1_, a2_] := 
  Integrate[jointf1[x1, x2], {x2, 0, a2}, {x1, 0, a1}, 
   Assumptions -> 
    0 < thB < 1 && 0 < xi < 1 && 0 < xbar < 1 && 0 < abarS < 1];
Plot3D[jointf1[x1, x2, 0.5, 0.91], {x1, 0, 1}, {x2, 0, 1}]

Thank you for your reply.

Let me explain: abarS, xbar and c are expressions of thB and xi. And jointf1 is defined as a function of four variables.

But when I try to plot the jointf1 with thB=0.5, xi=0.91, the code gives empty plot. Is there a problem in my code?

Thanks again!