Evaluation of integral with Hypergeometric1F1

Question

I solved the following inhomogeneous second order differential equation:

SolGeneralNonHomo = 
 DSolve[1/2*\[Sigma]X^2*V1''[x] + ((a0*\[Delta])/a1 - a1*x)*V1'[x] - 
     em*V1[x] + (2*x - 2*ep*(a0 + a1*x)* \[Lambda])/(em - ep) == 0, 
   V1[x], x] // FullSimplify

I want to plot one possible solution. So first I specify numerical values for the parameters:

\[Delta] = 0.99;
\[Theta] = 1.5;
\[Epsilon] = 1.5;
\[CurlyPhi] = 5.0;
\[Xi] = 0.02;
\[Lambda] = 0.5;
\[Sigma]X = 0.20;
\[Mu]A = 0.15;
\[Sigma]A = 0.20;
rn = \[Delta] + \[Mu]A - \[Sigma]A^2;
em = (\[Delta]*\[Theta] - 
   Sqrt[\[Delta]^2*\[Theta]^2 + 
    4*\[Theta]*(\[Epsilon] - 1)*(1 + \[CurlyPhi])])/(2*\[Theta]);
ep = (\[Delta]*\[Theta] + 
   Sqrt[\[Delta]^2*\[Theta]^2 + 
    4*\[Theta]*(\[Epsilon] - 1)*(1 + \[CurlyPhi])])/(2*\[Theta]);
a1 = (-\[Delta]*\[Theta] + 
   Sqrt[\[Delta]^2*\[Theta]^2 + 
    4*\[Theta]*(\[Epsilon] - 1)*(1 + \[CurlyPhi])])/(2*\[Theta]);
a0 = -((\[Xi] + rn - 
     1/2*\[Sigma]X^2 - \[Sigma]A*\[Sigma]X)/(\[Delta] + a1))*a1;

Then I assign the solution with these values to a new variable:

p[x_] = V1[x] /. SolGeneralNonHomo[[1]]

This yields:

9.9501 (0.528212 + 1.00669 x) C[1] + 
 C[2] Hypergeometric1F1[-0.5, 1/2, 
   24.7511 (0.528212 + 1.00669 x)^2] + 
 9.9501 (0.528212 + 1.00669 x) Inactive[Integrate][-((
    0.669635 Hypergeometric1F1[-0.5, 1/2, 
      24.7511 (-0.528212 - 1.00669 K[1])^2] (-0.531771 - 
       1.11022*10^-16 K[1]))/(
    0.201003 Hypergeometric1F1[-0.5, 1/2, 
       24.7511 (-0.528212 - 1.00669 K[1])^2] + 
     9.9501 Hypergeometric1F1[0.5, 3/2, 
       24.7511 (-0.528212 - 1.00669 K[1])^2] (0.528212 + 
        1.00669 K[1]) (0.528212 + 1.00669 K[1.]))), {K[1], 1, x}] + 
 Hypergeometric1F1[-0.5, 1/2, 
   24.7511 (0.528212 + 1.00669 x)^2] Inactive[Integrate][(
   6.66293 (-0.531771 - 1.11022*10^-16 K[2]) (0.528212 + 
      1.00669 K[2.]))/(
   0.201003 Hypergeometric1F1[-0.5, 1/2, 
      24.7511 (-0.528212 - 1.00669 K[2])^2] + 
    9.9501 Hypergeometric1F1[0.5, 3/2, 
      24.7511 (-0.528212 - 1.00669 K[2])^2] (0.528212 + 
       1.00669 K[2]) (0.528212 + 1.00669 K[2.])), {K[2], 1, x}]

Is there a reason why Mathematica can't evaluate this expression? For example, the following:

p[0.2]

Outputs this rather than a number:

7.2591 C[1] - 22873.1 C[2] + 
 7.2591 Inactive[Integrate][-((
    0.669635 Hypergeometric1F1[-0.5, 1/2, 
      24.7511 (-0.528212 - 1.00669 K[1])^2] (-0.531771 - 
       1.11022*10^-16 K[1]))/(
    0.201003 Hypergeometric1F1[-0.5, 1/2, 
       24.7511 (-0.528212 - 1.00669 K[1])^2] + 
     9.9501 Hypergeometric1F1[0.5, 3/2, 
       24.7511 (-0.528212 - 1.00669 K[1])^2] (0.528212 + 
        1.00669 K[1]) (0.528212 + 1.00669 K[1.]))), {K[1], 1, 0.2}] - 
 22873.1 Inactive[Integrate][(
   6.66293 (-0.531771 - 1.11022*10^-16 K[2]) (0.528212 + 
      1.00669 K[2.]))/(
   0.201003 Hypergeometric1F1[-0.5, 1/2, 
      24.7511 (-0.528212 - 1.00669 K[2])^2] + 
    9.9501 Hypergeometric1F1[0.5, 3/2, 
      24.7511 (-0.528212 - 1.00669 K[2])^2] (0.528212 + 
       1.00669 K[2]) (0.528212 + 1.00669 K[2.])), {K[2], 1, 0.2}]

Answer 1

$Version

(* "13.2.1 for Mac OS X x86 (64-bit) (January 27, 2023)" *)

Clear["Global`*"]

δ = 99/100;
θ = 3/2;
ϵ = 3/2;
φ = 5;
ξ = 1/50;
λ = 1/2;
σX = 1/5;
μA = 3/20;
σA = 1/5;
rn = δ + μA - σA^2;
em = (δ*θ - Sqrt[δ^2*θ^2 + 4*θ*(ϵ - 1)*(1 + φ)])/(2*θ);
ep = (δ*θ + Sqrt[δ^2*θ^2 + 4*θ*(ϵ - 1)*(1 + φ)])/(2*θ);
a1 = (-δ*θ + Sqrt[δ^2*θ^2 + 4*θ*(ϵ - 1)*(1 + φ)])/(2*θ);
a0 = -((ξ + rn - 1/2*σX^2 - σA*σX)/(δ + a1))*a1;

Solving,

SolGeneralNonHomo = 
 DSolve[1/2*σX^2*V1''[x] + ((a0*δ)/a1 - a1*x)*V1'[x] - 
      em*V1[x] + (2*x - 2*ep*(a0 + a1*x)*λ)/(em - ep) == 0, V1[x], 
    x][[1]] // FullSimplify

(* {V1[x] -> (2 (5300 (99 + Sqrt[89801]) + 
        Sqrt[89801 (99 + Sqrt[89801])] (5247 + 10000 x) C[1] + 
        50 (89801 + 99 Sqrt[89801]) E^((5247 + 10000 x)^2/(
         10000 (99 + Sqrt[89801]))) C[2]) - 
     89801^(1/4) Sqrt[(89801 + 99 Sqrt[89801]) π] (5247 + 10000 x) C[
      2] Erfi[(Sqrt[1/2 (-99 + Sqrt[89801])] (5247 + 10000 x))/
       20000])/(100 (89801 + 99 Sqrt[89801]))} *)

There are two arbitrary constants (C[1] and C[2]). Assign initial conditions to resolve these. For example, if V1[0] == 1, V1'[0] == 1 then

SolGeneralNonHomo2 = 
 DSolve[{1/2*σX^2*V1''[x] + ((a0*δ)/a1 - a1*x)*V1'[x] - 
      em*V1[x] + (2*x - 2*ep*(a0 + a1*x)*λ)/(em - ep) == 0, 
    V1[0] == 1, V1'[0] == 1}, V1[x], x] // FullSimplify

(* {{V1[x] -> (E^(-(27531009/(
       10000 (99 + Sqrt[89801])))) (100 89801^(
         1/4) (-52933468853 + 321884153 Sqrt[89801]) E^((5247 + 10000 x)^2/(
         10000 (99 + Sqrt[89801]))) + 
        100 89801^(1/4) E^(27531009/(
         10000 (99 + Sqrt[89801]))) (53 (2676159601 + 10870299 Sqrt[89801]) + 
           898010000 (99 + Sqrt[89801]) x) - (-426824153 + 
           1060000 Sqrt[89801]) Sqrt[(89801 + 99 Sqrt[89801]) π] (5247 + 
           10000 x) Erfi[(5247 Sqrt[1/2 (-99 + Sqrt[89801])])/
          20000] + (-426824153 + 
           1060000 Sqrt[89801]) Sqrt[(89801 + 99 Sqrt[89801]) π] (5247 + 
           10000 x) Erfi[(Sqrt[1/2 (-99 + Sqrt[89801])] (5247 + 10000 x))/
          20000]))/(89801000000 89801^(1/4) (99 + Sqrt[89801]))}} *)

p[x_] = V1[x] /. SolGeneralNonHomo2[[1]];

p[0.2]

(* -1.52323 *)

Plot[p[x], {x, 0, 1/4}]