Numerical Double Integral with symbolic bounds to normalize function

Question

I'm trying to use mathematica to normalize this trial wave function which requires finding the following integral and solving for A so that the total propability equals 1. 'a' is the constant width of the potential well and alpha is the variational parameter. I'm having a lot of dificulty with mathematica accepting non-numeric bounds such as a and while still managing to get numerical results but it also can't be done purely symbolically without encountering a ton of erf functions.

8*A^2*Integrate[((1 - x^2*y^2)/a^4/E^(\[Alpha]*((x^2 + y^2)/a^2)))^2, {y, 0, a}, {x, y, a}] + 8*A^2*Integrate[((1 - x^2)/a^4/E^(\[Alpha]*((x^2 + y^2)/a^2)))^2, {y, 0, a}, {x, a, Infinity}]

f1[(x_)?NumericQ, (y_)?NumericQ] := NIntegrate[((1 - x^2*y^2)/a^4 E^(\[Alpha]*((x^2 + y^2)/a^2)))^2, {x, y, a}]

NIntegrate[f1[y], {y, 0, a}]

NIntegrate::nlim: y = a is not a valid limit of integration.


N[8*A^2*Integrate[((1 - x^2*y^2)/a^4/
   E^(\[Alpha]*((x^2 + y^2)/a^2)))^2, {y, 0, a}, 
 {x, y, a}] + 8*A^2*Integrate[
 ((1 - x^2)/a^4/E^(\[Alpha]*((x^2 + y^2)/a^2)))^
  2, {y, 0, a}, {x, a, Infinity}]]

((1/(a^6*\[Alpha]^5))*0.001953125*A^2*
    (8.*a^4*\[Alpha]*(-32.*\[Alpha]^2 + 
            a^4*(3. + 4.*\[Alpha])^2) - 10.026513098524001*
         2.718281828459045^(2.*\[Alpha])*a^4*Sqrt[\[Alpha]]*
         (-32.*\[Alpha]^2 + 3.*a^4*(3. + 4.*\[Alpha]))*
         Erf[1.4142135623730951*Sqrt[\[Alpha]]] + 
       3.141592653589793*2.718281828459045^
           (4.*\[Alpha])*(9.*a^8 - 32.*a^4*\[Alpha]^2 + 
            256.*\[Alpha]^4)*Erf[1.4142135623730951*
               Sqrt[\[Alpha]]]^2))/2.718281828459045^
   (4.*\[Alpha]) + ((1/(a^6*\[Alpha]^3))*0.0625*A^2*
    Erf[1.4142135623730951*Sqrt[\[Alpha]]]*
    (5.0132565492620005*a^2*Sqrt[\[Alpha]]*
         (-8.*\[Alpha] + a^2*(3. + 4.*\[Alpha])) + 
       3.141592653589793*2.718281828459045^
           (2.*\[Alpha])*(3.*a^4 - 8.*a^2*\[Alpha] + 
     16.*\[Alpha]^2)*
         Erfc[1.4142135623730951*Sqrt[\[Alpha]]]))/
 2.718281828459045^(2.*\[Alpha])

Answer 1

Analytical integration yields good results when applying assumptions.

int = 8*A^2*
        Integrate[((1 - x^2*y^2)/a^4/E^(\[Alpha]*((x^2 + y^2)/a^2)))^2, {y,
  0, a}, {x, y, a}, Assumptions -> a > 0 && \[Alpha] > 0] + 
     8*A^2*Integrate[((1 - x^2)/a^4/
   E^(\[Alpha]*((x^2 + y^2)/a^2)))^2, {y, 0, a}, {x, a, Infinity},
 Assumptions -> a > 0 && \[Alpha] > 0]

Result is a long expression. Solve for A:

AA[a_, \[Alpha]_] = A /. Solve[1 == int, A]

Get A for a given a = 1:

AA[1, \[Alpha]] // FullSimplify[# \[Alpha] > 0] &

Since A^2 is used, you have two solutions:

Plot3D[Evaluate@AA[a, \[Alpha]], {a, 0, 5}, {\[Alpha], 0, 4}, 
    PlotStyle -> {Blue, Red}]