h=0.7; cspeed = 299792.458; Omega[CapitalLambda]0 = 0.7; Omegam0 = 0.3; Omega0 = 1.0002;
Yp = Interval[0.2534 + 0.0083 {-1, 1}]; alpha = Interval[0.17 + 0.03 {-1, 1}]; logvc0 =  Interval[1.58 + 0.05 {-1, 1}]; logvc1 = Interval[3.14 + 0.38 {-1, 1}]; beta = Interval[-0.50 + 0.18 {-1, 1}];
var1 = Interval[10.79 + 0.01 {-1, 1}]; phi1 = Interval[-3.31 + 0.20 {-1, 1}]; phi2 = Interval[-2.01 + 0.28 {-1, 1}]; alpha1 = Interval[-1.69 + 0.10 {-1, 1}]; alpha2 = Interval[-0.79 + 0.04 {-1, 1}];
[Alpha]0 = Interval[-1.412 + {-0.105, 0.020}]; [Alpha]a = Interval[0.731 + {-0.296, 0.344}];
[Delta]0 = Interval[3.508 + {-0.369, 0.087}]; [Delta]a = Interval[2.608 + {-1.261, 2.446}]; [Delta]z = Interval[-0.043 + {-0.071, 0.958}];
[Gamma]0 = Interval[0.316 + {-0.012, 0.076}]; [Gamma]a =  Interval[1.319 + {-1.261, 0.584}]; [Gamma]z = Interval[0.279 + {-0.081, 0.256}];
[Epsilon]0 = Interval[-1.777 + {-0.146, 0.133}]; [Epsilon]a = Interval[-0.006 + {-0.361, 0.113}]; [Epsilon]z = Interval[0.0 + {-0.104, 0.003}]; [Epsilon]a2 = Interval[-0.119 + {-0.012, 0.061}];
M10 = Interval[11.514 + {-0.009, 0.053}]; M1a = Interval[-1.793 + {-0.330, 0.315}]; M1z = Interval[-0.251 + {-0.125, 0.012}];

[Alpha][z_?NumericQ] := [Alpha][z] = 1./(1. + z);
distance[z_?NumericQ] := distance[z] = NIntegrate[(Omegam0*(1 + u)^3 + Omega[CapitalLambda]0 + (Omega0 - Omegam0 - Omega[CapitalLambda]0)*(1 + u)^2)^(-1/2), {u, 0, z}];
F[z_?NumericQ] := F[z] = (Sin[distance[z]*(Omega0 - 1.)^(1/2)]/(distance[z]*(Omega0 - 1.)^(1/2)));
ComovingVolume[z_?NumericQ] := ComovingVolume[z] = (4*[Pi])/3*(cspeed/100)^3*(distance[z])^3*((3*(1 - (Sin[2*distance[z]*(Omega0 - 1.)^(1/2)]/(2*distance[z]*(Omega0 - 1.)^(1/2)))))/(2*(Omega0 - 1)*(distance[z])^2));
dComovingVolumedz[z_?NumericQ] := dComovingVolumedz[z] = Derivative[1][ComovingVolume][z] // N;
[Nu]func[z_?NumericQ] := [Nu]func[z] = Exp[-4 *([Alpha][z] )^2];
[Alpha]func[z_?NumericQ, [Alpha]0_, [Alpha]a_] := [Alpha]func[z, [Alpha]0, [Alpha]a] = [Alpha]0 + ([Alpha]a *([Alpha][z] - 1))*[Nu]func[z];
 [Delta]func[z_?NumericQ, [Delta]0_, [Delta]a_, [Delta]z_] := [Delta]func[z, [Delta]0, [Delta]a, [Delta]z] = [Delta]0 + ([Delta]a *(
[Alpha][z] - 1) + [Delta]z* z) *[Nu]func[z];
[Gamma]func[z_?NumericQ, [Gamma]0_, [Gamma]a_, [Gamma]z_] := [Gamma]func[z, [Gamma]0, [Gamma]a, [Gamma]z] = [Gamma]0 + ([Gamma]a (
[Alpha][z] - 1) + [Gamma]z *z) *[Nu]func[z];
log10[Epsilon][z_?NumericQ, [Epsilon]0_, [Epsilon]a_, [Epsilon]z_,
[Epsilon]a2_] := log10[Epsilon][z, [Epsilon]0, [Epsilon]a, [Epsilon]z, [Epsilon]a2] = [Epsilon]0 + ([Epsilon]a*([Alpha][z] - 1) + [Epsilon]z*z)*[Nu]func[z] + [Epsilon]a2*([Alpha][z] - 1);
log10M1[z_?NumericQ, M10_, M1a_, M1z_] := log10M1[z, M10, M1a, M1z] =
   Log10[h] + M10 + (M1a*([Alpha][z] - 1) + M1z*z)*[Nu]func[z];
fun[y_?NumericQ, z_?NumericQ, [Alpha]0_, [Alpha]a_, [Delta]0_, [Delta]a_, [Delta]z_, [Gamma]0_, [Gamma]a_, [Gamma]z_] :=
  fun[y, z, [Alpha]0, [Alpha]a, [Delta]0, [Delta]a, [Delta]z, [Gamma]0, [Gamma]a, [Gamma]z] = -Log10[10^( ([Alpha]func[z, [Alpha]0, [Alpha]a])*y) + 1] + [Delta]func[z, [Delta]0, [Delta]a, [Delta]z] *(Log10[1 + Exp[y]])^ [Gamma]func[z, [Gamma]0, [Gamma]a, [Gamma]z]/(1 + Exp[10^-y]);
MStar[z_?NumericQ, M_?NumericQ, [Epsilon]0_, [Epsilon]a_, [Epsilon]z_, [Epsilon]a2_, M10_, M1a_, M1z_, [Alpha]0_, [Alpha]a_, [Delta]0_, [Delta]a_, [Delta]z_, [Gamma]0_, [Gamma]a_, [Gamma]z_] :=
  MStar[z, M, [Epsilon]0, [Epsilon]a, [Epsilon]z, [Epsilon]a2, M10, M1a, M1z, [Alpha]0, [Alpha]a, [Delta]0, [Delta]a, [Delta]z, [Gamma]0, [Gamma]a, [Gamma]z] = log10[Epsilon][z, [Epsilon]0, [Epsilon]a, [Epsilon]z, [Epsilon]a2] + log10M1[z, M10, M1a, M1z] + fun[M - log10M1[z, M10, M1a, M1z], z, [Alpha]0, [Alpha]a, [Delta]0, [Delta]a, [Delta]z, [Gamma]0, [Gamma]a, [Gamma]z] - fun[0, z, [Alpha]0, [Alpha]a, [Delta]0, [Delta]a, [Delta]z, [Gamma]0, [Gamma]a, [Gamma]z];
dLogMStardLogM[z_?NumericQ, M_?NumericQ, [Epsilon]0_, [Epsilon]a_, [Epsilon]z_, [Epsilon]a2_, M10_, M1a_, M1z_, [Alpha]0_, [Alpha]a_, [Delta]0_, [Delta]a_, [Delta]z_, [Gamma]0_, [Gamma]a_, [Gamma]z_] := dLogMStardLogM[z, M, [Epsilon]0, [Epsilon]a, [Epsilon]z, [Epsilon]a2, M10, M1a, M1z, [Alpha]0, [Alpha]a, [Delta]0, [Delta]a, [Delta]z, [Gamma]0, [Gamma]a, [Gamma]z] = Derivative[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0][MStar][z, M, [Epsilon]0, [Epsilon]a, [Epsilon]z, [Epsilon]a2, M10, M1a, M1z, [Alpha]0, [Alpha]a, [Delta]0, [Delta]a, [Delta]z, [Gamma]0, [Gamma]a, [Gamma]z];
LocalSMFAllGalaxiesWSB[MS_?NumericQ, var1_, phi1_, phi2_, alpha1_, alpha2_] := LocalSMFAllGalaxiesWSB[MS, var1, phi1, phi2, alpha1, alpha2] =
   Log[10.]*Exp[-(10.^(MS - var1 - Log10[h]))]*(10.^phi1*(10.^(MS - var1 - Log10[h]))^(alpha1 + 1.) + 10.^phi2*(10.^(MS - var1 - Log10[h]))^(alpha2 + 1.));
LocalSMF[z_?NumericQ, M_?NumericQ, var1_, phi1_, phi2_, alpha1_, alpha2_, [Epsilon]0_, [Epsilon]a_, [Epsilon]z_, [Epsilon]a2_, M10_, M1a_, M1z_, [Alpha]0_, [Alpha]a_, [Delta]0_, [Delta]a_, [Delta]z_, [Gamma]0_, [Gamma]a_, [Gamma]z_] := LocalSMF[z, M, var1, phi1, phi2, alpha1, alpha2, [Epsilon]0, [Epsilon]a, [Epsilon]z, [Epsilon]a2, M10, M1a, M1z, [Alpha]0, [Alpha]a, [Delta]0, [Delta]a, [Delta]z, [Gamma]0, [Gamma]a, [Gamma]z] = LocalSMFAllGalaxiesWSB[MStar[z, M, [Epsilon]0, [Epsilon]a, [Epsilon]z, [Epsilon]a2, M10, M1a, M1z, [Alpha]0, [Alpha]a, [Delta]0, [Delta]a, [Delta]z, [Gamma]0, [Gamma]a, [Gamma]z], var1, phi1, phi2, alpha1, alpha2];
dNgalaxiesdzWSB[z_, M_, var1_, phi1_, phi2_, alpha1_, alpha2_, [Epsilon]0_, [Epsilon]a_, [Epsilon]z_, [Epsilon]a2_, M10_, M1a_, M1z_, [Alpha]0_, [Alpha]a_, [Delta]0_, [Delta]a_, [Delta]z_, [Gamma]0_, [Gamma]a_, [Gamma]z_] := dNgalaxiesdzWSB[z, M, var1, phi1, phi2, alpha1, alpha2, [Epsilon]0, [Epsilon]a, [Epsilon]z, [Epsilon]a2, M10, M1a, M1z, [Alpha]0, [Alpha]a, [Delta]0, [Delta]a, [Delta]z, [Gamma]0, [Gamma]a, [Gamma]z] = 7748./(4.*[Pi])*([Pi]/180.)^2.*dComovingVolumedz[z]*Integrate[LocalSMF[z, p, var1, phi1, phi2, alpha1, alpha2, [Epsilon]0, [Epsilon]a, [Epsilon]z, [Epsilon]a2, M10, M1a, M1z, [Alpha]0, [Alpha]a, [Delta]0, [Delta]a, [Delta]z, [Gamma]0, [Gamma]a, [Gamma]z], {p, M, 15.}];

LogPlot[{
Min[dNgalaxiesdzWSB[0.03, M, var1, phi1, phi2, alpha1, alpha2, [Epsilon]0, [Epsilon]a, [Epsilon]z, [Epsilon]a2, M10, M1a, M1z, [Alpha]0, [Alpha]a, [Delta]0, [Delta]a, [Delta]z, [Gamma]0, [Gamma]a, [Gamma]z]],
Max[dNgalaxiesdzWSB[0.03, M, var1, phi1, phi2, alpha1, alpha2, [Epsilon]0, [Epsilon]a, [Epsilon]z, [Epsilon]a2, M10, M1a, M1z, [Alpha]0, [Alpha]a, [Delta]0, [Delta]a, [Delta]z, [Gamma]0, [Gamma]a, [Gamma]z]]},
{M, 10.25, 12.95},
PlotRange -> {10^-3, 10^7.4},
 Filling -> {1 -> {{2}, LightGreen}},
 PlotStyle -> {{Green, Thickness[0.005]}, {Green, Thickness[0.005]}},
 Frame -> True, FrameLabel -> {Style["X", FontSize -> 26], Style["Y", FontSize -> 26]},
 FrameTicksStyle -> Directive[FontSize -> 26]]