Hi All, I have a question which has bothered me for a whole week. I defined some global functions and then calculated two specific function values based on them. I can plot the first one before running the calculation for the second function, but not after I running the calculation for the second function. To make my confusion clear, I copied my code below:

(*Define public parameters*)
c = 2.99792458*10^8;
kB = 1.38*10^-23;
\[Mu]0 = 4 \[Pi]*10^-7;
n1 = 1.4469;
n2 = 1.0;
a = 200*10^-9;
\[Omega][\[Lambda]_] := 2*\[Pi]*c/\[Lambda];
k0[\[Lambda]_] := 2*\[Pi]/\[Lambda];
JD[l_, x_] := BesselJ[l - 1, x] - l*BesselJ[l, x]/x;
KD[l_, x_] := -(BesselK[l - 1, x] + BesselK[l + 1, x])/2;
s11[\[Lambda]_] := (1/(h11[\[Lambda]]*a)^2 + 
     1/(q11[\[Lambda]]*a)^2)/(JD[1, h11[\[Lambda]]*a]/h11[\[Lambda]]/
       a/BesselJ[1, h11[\[Lambda]]*a] + 
     KD[1, q11[\[Lambda]]*a]/q11[\[Lambda]]/a/
      BesselK[1, q11[\[Lambda]]*a]);
f[\[Lambda]_] := (1 + s11[\[Lambda]])^2/(1 - s11[\[Lambda]])^2;
fp[\[Lambda]_] := 2*(1 + s11[\[Lambda]])/(1 - s11[\[Lambda]]);
r[x_, y_] := Sqrt[x^2 + y^2];
e = 1.602176*10^-19;
me = 9.109381*10^-31;
\[Omega]\[Alpha] = {2 \[Pi]*384.2304844685*10^12, 
   2 \[Pi]*377.107463380*10^12};(*D2 and D1 lines*)
f\[Alpha] = {0.69577, 0.34231};(*Oscillator damping coefficients*)
\[Gamma]\[Alpha] = {38.117*10^6, 36.129*10^6};(*natual line widths*)
\[Alpha]1[\[Lambda]_] := 
  e^2/me*Sum[
    f\[Alpha][[
      i]]*(\[Omega]\[Alpha][[
         i]]^2 - \[Omega][\[Lambda]]^2)/((\[Omega]\[Alpha][[
           i]]^2 - \[Omega][\[Lambda]]^2)^2 + \[Gamma]\[Alpha][[
          i]]^2*\[Omega][\[Lambda]]^2), {i, 1, 2} ];
(*Blue-detuned parameters*)
P1 = 0.029;
\[Lambda]1 = 700*10^-9;
\[Beta]11 = 10.1477*10^6;(*propagation constant*)
q11[\[Lambda]_] := Sqrt[\[Beta]11^2 - k0[\[Lambda]]^2*n2^2];
h11[\[Lambda]_] := Sqrt[k0[\[Lambda]]^2*n1^2 - \[Beta]11^2];
Din1[\[Lambda]_] := (1 - 
      s11[\[Lambda]])*(1 + (1 - s11[\[Lambda]])*\[Beta]11^2/
        h11[\[Lambda]]^2)*((BesselJ[0, 
        h11[\[Lambda]]*a])^2 + (BesselJ[1, 
        h11[\[Lambda]]*a])^2) + (1 + 
      s11[\[Lambda]])*(1 + (1 + s11[\[Lambda]])*\[Beta]11^2/
        h11[\[Lambda]]^2)*((BesselJ[2, h11[\[Lambda]]*a])^2 - 
      BesselJ[1, h11[\[Lambda]]*a]*BesselJ[3, h11[\[Lambda]]*a]);
Dout1[\[Lambda]_] := (BesselJ[1, 
     h11[\[Lambda]]*
      a])^2*((1 - 
         s11[\[Lambda]])*(1 - (1 - s11[\[Lambda]])*\[Beta]11^2/
           q11[\[Lambda]]^2)*((BesselK[0, 
           q11[\[Lambda]]*a])^2 - (BesselK[1, 
           q11[\[Lambda]]*a])^2) + (1 + 
         s11[\[Lambda]])*(1 - (1 + s11[\[Lambda]])*\[Beta]11^2/
           q11[\[Lambda]]^2)*((BesselK[2, q11[\[Lambda]]*a])^2 - 
         BesselK[1, q11[\[Lambda]]*a]*
          BesselK[3, q11[\[Lambda]]*a]))/(BesselK[1, 
      q11[\[Lambda]]*a])^2;
A1[\[Lambda]_] := 
  Sqrt[4*\[Mu]0*\[Omega][\[Lambda]]*
     P1/\[Pi]/a^2/\[Beta]11]*(Din1[\[Lambda]] + Dout1[\[Lambda]])^(-1/
    2);
Alin1[\[Lambda]_] := \[Sqrt]2*A1[\[Lambda]];
u1[\[Lambda]_] := 
  2*h11[\[Lambda]]^2/\[Beta]11^2/(1 - s11[\[Lambda]])^2;
w1[\[Lambda]_] := 
  2*q11[\[Lambda]]^2/\[Beta]11^2/(1 - s11[\[Lambda]])^2;
gin1[\[Lambda]_] := (Abs[Alin1[\[Lambda]]])^2/2/u1[\[Lambda]];
gout1[\[Lambda]_] := (Abs[
      Alin1[\[Lambda]]])^2*(BesselJ[1, h11[\[Lambda]]*a])^2/2/
     w1[\[Lambda]]/(BesselK[1, q11[\[Lambda]]*a])^2;
(*circularly-polarized*)
I1in[r_] := 
  2 gin1[\[Lambda]1]*((BesselJ[0, h11[\[Lambda]1]*r])^2 + 
     u1[\[Lambda]1]*(BesselJ[1, h11[\[Lambda]1]*r])^2 + 
     f[\[Lambda]1]*(BesselJ[2, h11[\[Lambda]1]*r])^2);
I1out[r_] := 
  2 gout1[\[Lambda]1]*((BesselK[0, q11[\[Lambda]1]*r])^2 + 
     w1[\[Lambda]1]*(BesselK[1, q11[\[Lambda]1]*r])^2 + 
     f[\[Lambda]1]*(BesselK[2, q11[\[Lambda]1]*r])^2);
F2[x_, y_] := (I1in[r[x, y]]/I1out[a]) /; Sqrt[x^2 + y^2] <= a
F2[x_, y_] := (I1out[r[x, y]]/I1out[a]) /; Sqrt[x^2 + y^2] > a
Fig6 = Plot3D[F2[x, y], {x, -4 a, 4 a}, {y, -4 a, 4 a}, 
  PlotRange -> All]
U1[x_, y_] := 
  1000*(-\[Alpha]1[\[Lambda]1]*F2[x, y]/4)*2/
    kB;(*potential energy=kB*T/2,unit in mK*)
Fig7 = Plot3D[U1[x, y], {x, a, 4*a}, {y, a, 4*a}, 
  PlotRange -> Automatic]
(*Red-detuned parameters*)
P2 = 0.03;
\[Lambda]2 = 1060*10^-9;
\[Beta]11 = 6.04866*10^6;(*propagation constant*)
q11[\[Lambda]_] := Sqrt[\[Beta]11^2 - k0[\[Lambda]]^2*n2^2];
h11[\[Lambda]_] := Sqrt[k0[\[Lambda]]^2*n1^2 - \[Beta]11^2];
Din2[\[Lambda]_] := (1 - 
      s11[\[Lambda]])*(1 + (1 - s11[\[Lambda]])*\[Beta]11^2/
        h11[\[Lambda]]^2)*((BesselJ[0, 
        h11[\[Lambda]]*a])^2 + (BesselJ[1, 
        h11[\[Lambda]]*a])^2) + (1 + 
      s11[\[Lambda]])*(1 + (1 + s11[\[Lambda]])*\[Beta]11^2/
        h11[\[Lambda]]^2)*((BesselJ[2, h11[\[Lambda]]*a])^2 - 
      BesselJ[1, h11[\[Lambda]]*a]*BesselJ[3, h11[\[Lambda]]*a]);
Dout2[\[Lambda]_] := (BesselJ[1, 
     h11[\[Lambda]]*
      a])^2*((1 - 
         s11[\[Lambda]])*(1 - (1 - s11[\[Lambda]])*\[Beta]11^2/
           q11[\[Lambda]]^2)*((BesselK[0, 
           q11[\[Lambda]]*a])^2 - (BesselK[1, 
           q11[\[Lambda]]*a])^2) + (1 + 
         s11[\[Lambda]])*(1 - (1 + s11[\[Lambda]])*\[Beta]11^2/
           q11[\[Lambda]]^2)*((BesselK[2, q11[\[Lambda]]*a])^2 - 
         BesselK[1, q11[\[Lambda]]*a]*
          BesselK[3, q11[\[Lambda]]*a]))/(BesselK[1, 
      q11[\[Lambda]]*a])^2;
A2[\[Lambda]_] := 
  Sqrt[4*\[Mu]0*\[Omega][\[Lambda]]*
     P2/\[Pi]/a^2/\[Beta]11]*(Din2[\[Lambda]] + Dout2[\[Lambda]])^(-1/
    2);
Alin2[\[Lambda]_] := \[Sqrt]2*A2[\[Lambda]];
u2[\[Lambda]_] := 
  2*h11[\[Lambda]]^2/\[Beta]11^2/(1 - s11[\[Lambda]])^2;
w2[\[Lambda]_] := 
  2*q11[\[Lambda]]^2/\[Beta]11^2/(1 - s11[\[Lambda]])^2;
gin2[\[Lambda]_] := (Abs[Alin2[\[Lambda]]])^2/2/u2[\[Lambda]];
gout2[\[Lambda]_] := (Abs[
      Alin2[\[Lambda]]])^2*(BesselJ[1, h11[\[Lambda]]*a])^2/2/
     w2[\[Lambda]]/(BesselK[1, q11[\[Lambda]]*a])^2;
I3in[r_] := 
  2 gin2[\[Lambda]2]*((BesselJ[0, h11[\[Lambda]2]*r])^2 + 
     u2[\[Lambda]2]*(BesselJ[1, h11[\[Lambda]2]*r])^2 + 
     f[\[Lambda]2]*(BesselJ[2, h11[\[Lambda]2]*r])^2);
I3out[r_] := 
  2 gout2[\[Lambda]2]*((BesselK[0, q11[\[Lambda]2]*r])^2 + 
     w2[\[Lambda]2]*(BesselK[1, q11[\[Lambda]2]*r])^2 + 
     f[\[Lambda]2]*(BesselK[2, q11[\[Lambda]2]*r])^2);
F3[x_, y_] := (I3in[r[x, y]]/I3out[a]) /; Sqrt[x^2 + y^2] <= a
F3[x_, y_] := (I3out[r[x, y]]/I3out[a]) /; Sqrt[x^2 + y^2] > a
Fig8 = Plot3D[F3[x, y], {x, -4 a, 4 a}, {y, -4 a, 4 a}, 
  PlotRange -> All]
U2[x_, y_] := 
  1000*(-\[Alpha]1[\[Lambda]2]*F3[x, y]/4)*2/
    kB;(*potential energy=kB*T/2,unit in mK*)
Fig9 = Plot3D[U2[x, y], {x, a, 4*a}, {y, a, 4*a}, 
  PlotRange -> Automatic]

I am not sure whether you need to copy the code into Mathematica to see it clearly. The problem is I can only plot out U1[x,y] before running the calculation for U2[x,y] but not after. I also attached my .nb file with this posting. Thank you very much for your help and look forward to any advice. Best, Dienter code here

Attachments: