How to solve the following complex integral?

Question

I am currently trying to reproduce a semi-analytical solution, where $U_0$ satisfies the PDE under the Black-Scholes model, $\mathcal{M}$ represents some differential operators, such as derivatives with respect to time $t$. Given such a complex recursive form involving the integrals of $U_n$, can Mathematica handle it? I found myself quite foolish as I directly used Integrate and failed to obtain a solution.

My code is:

T = 1;
r = 0.04;
kappa = 1;
theta = 0.25;
sigma = 0.1;

eta[\[Tau]_, x_, v_] := 
 Exp[-1/8 v (2 r/v + 1)^2 \[Tau] + 1/2 (2 r/v + 1) x]
dmPlus[t_, x_, v_] := (-x + (r + v/2) (T - t))/Sqrt[v (T - t)]
dmMinus[t_, x_, v_] := (-x + (r - v/2) (T - t))/Sqrt[v (T - t)]
dmR[t_, x_, v_] := dmPlus[t, x, v] - 2 r/Sqrt[v] Sqrt[T - t]
G[t_, x_, \[Xi]_, v_] := 
 1/(Sqrt[2 \[Pi] v t])*(Exp[-(x - \[Xi])^2/(2 v t)] + 
    Exp[-(x + \[Xi])^2/(2 v t)] + (2 r/v + 1) Integrate[
      Exp[-(x + \[Xi] + \[Eta])^2/(2 v t) + 
        0.5 (2 r/v + 1) \[Eta]], {\[Eta], 0, \[Infinity]}, 
      Assumptions -> 
       1/(t v) > 0 && (0.04` t - 1.` x - 1.` \[Xi])/(t v) < -0.5])
normal[z_] := CDF[NormalDistribution[0, 1], z]

U0[t_, x_, v_] := 
 Exp[x] Exp[-r (T - t)] normal[-dmMinus[t, x, v]] - 
  normal[-dmPlus[t, x, v]] + 
  v/(2 r) (Exp[-r (T - t)] normal[dmPlus[t, x, v]] - 
     Exp[x + 2 x r/v] normal[dmR[t, x, v]])
MU0[t_, x_, v_] := 
 kappa (theta - v) D[U0[t, x, v], v] + 
  0.5 sigma^2 v D[D[U0[t, x, v], v], v]
Uhat[\[Tau]_, x_, v_] := 
 Integrate[
  Integrate[
   eta[s, \[Xi], v] *MU0[T - s, \[Xi], v] *
    G[\[Tau] - s, x, \[Xi], v], {\[Xi], 0, \[Infinity]}, 
   Assumptions -> \[Tau] > s > 0 && x > 0 && v > 0], {s, 0, \[Tau]}]