Some necessary constants :

e = 1.602176487 10^(-19)
m0 = 9.10938215 10^(-31)
e0 = 8.85418781762 10^(-12)
hbar = 1.054571628 10^(-34)

Bohr radius :

a0 = ((4 Pi e0 (hbar^2))/(m0 (e^2)))

The 2pz function in spherical coordinates:

ψ2pz[r_, θ_, ϕ_] := (1/(4 (Sqrt[2 Pi])))((1 )/(a0^(3/2)))(r/a0) Exp[-(r/(2 a0))]Cos[θ]

Is this function normalized?

Integrate[ψ2pz[r,θ,ϕ] ψ2pz[r,θ,ϕ](r^2)Sin[θ],{r, 0, Infinity},{θ, 0, Pi},{ϕ, 0, 2 Pi}]

The 2pz function in cartesian coordinates :

ψ2pzCc[x_, y_, z_] := (1/(4 (Sqrt[2 Pi])))((1 )/(a0^(3/2)))((Sqrt[(x^2) + (y^(2 ))+ (z^2)])/a0) Exp[-((Sqrt[(x^2) + (y^(2 ))+ (z^2)])/(2 a0))](z/(Sqrt[(x^2) + (y^(2 ))+ (z^2)]))

The 2s function in spherical coordinates:

ψ2s[r_, θ_, ϕ_] := (1/( 4 (Sqrt[ 2 Pi])))((1 )/(a0^(3/2)))Exp[-(r/( 2 a0))](2 - (r/a0))

Is this function normalized?

Integrate[ψ2s[r,θ,ϕ] ψ2s[r,θ,ϕ](r^2)Sin[θ],{r, 0, Infinity},{θ, 0, Pi},{ϕ, 0, 2 Pi}]

The 2s function in cartesian coordinates :

ψ2sCc[x_, y_, z_] := (1/( 4 (Sqrt[ 2 Pi])))((1 )/(a0^(3/2)))Exp[-((Sqrt[(x^2) + (y^(2 ))+ (z^2)])/( 2 a0))](2 - ((Sqrt[(x^2) + (y^(2 ))+ (z^2)])/a0))

Plot the 2s and the 2pz function, as well as the hybrid function as a function of z at x=y=0:

Plot[{ψ2sCc[0,0,z a0],ψ2pzCc[0,0,z a0],(1/(Sqrt[2]))(ψ2sCc[0,0,z a0]+ψ2pzCc[0,0,z a0])},{z,-7.5, 7.5}, PlotRange -> {{-7.5 , 7.5 },{-2×10^(14), 6×10^(14)}}]

Make a contour plot of the hybrid function :

ContourPlot[(1/(Sqrt[2]))(ψ2sCc[x a0,0,z a0]+ψ2pzCc[x a0,0,z a0]), {x, -10, 10}, {z, -10, 10},PlotRange -> {{-10, 10},{-10, 10},{-4 ×10^(14), 4 10^(14)}},ColorFunction->"Rainbow",Contours -> 30]

Make a contour plot of the 2pz function for comparison:

ContourPlot[ψ2pzCc[x a0,0,z a0], {x, -10, 10}, {z, -10, 10},PlotRange -> {{-10, 10},{-10, 10},{-4 ×10^(14), 4 10^(14)}},ColorFunction->"Rainbow",Contours -> 30]

ContourPlot[(1/2)(((ψ2sCc[x a0,0,z a0]+ψ2pzCc[x a0,0,z a0]))^2), {x, -10, 10}, {z, -10, 10},PlotRange -> {{-10, 10},{-10, 10},{0, (((3 10^(14)))^2)}},ColorFunction->"Rainbow",Contours -> 30]