(* Strategy: *)
(* 1. Make a grid of u and x values. *)
(* 2. For each x value, perform the luminosity integral over the u values. *)
(* 3. Plot the result vs. x. *)

(SubStar[τ])=1;
R0=1.5;
q=0;
β=1;

ux[x_,β_]:=Min[1/R0,1-((Abs[x])^(1/β))]

t[p_,z_]:=Piecewise[{{∞,Abs[p]≤1&&z≤(Sqrt[1-(p^2)])}},NIntegrate[((SubStar[τ])/(((p^2)+(Z^2)) (((1-(1/(Sqrt[(p^2)+(Z^2)]))))^β))),{Z,z,∞}]]

τ[μ_,r_]:=t[Re[(Sqrt[1-(μ^2)])] r,μ r]

f[r_]:=(r^(-q))

lum[x_,β_]:= NIntegrate[((f[1/u])/(((1-u))^(3 β))) Exp[-τ[-x/(((1-u))^β),1/u]],{u,0,ux[x,β]}]

(* u values for the grid *)
uArray[x_,β_]:=Table[Round[u,.01],{u,0,ux[x,β],.01}]

(* x values for the grid *)
xArray=Table[Round[x,0.01],{x,-1,1,0.01}];
xLum=xArray; (* initialize the luminosity array *)

(* loop over x values *)
For[j=1,j<201,j++,
x=xArray[[j]];
uarr=uArray[x,β]; (* initialize the x-specific u values *)
arr=uarr; (* initialize an array for values of the integrand *)
arr[[1]]=1;
(* loop over u values *)
For[i=2,i<ux[x,β]/0.01+1,i++,
temp=uarr[[i]];
(* calculate the integrand at each u value *)
arr[[i]]=((f[1/temp])/(((1-temp))^(3 β))) Exp[-τ[-x/(((1-temp))^β),1/temp]]; 
];
(* interpolate the integrand *)
func=Interpolation[Transpose[Join[{uarr},{arr}]]];
(* integrate the integrand for a certain value of x 
and store to the luminosity array*)
xLum[[j]]=NIntegrate[func[u],{u,0,ux[x,β]}];
]
xLum[[201]]=0;
(* interpolate the luminosity array *)
profile=Interpolation[Transpose[Join[{xArray},{xLum}]]]

Plot[profile[x],{x,-1,1}]