I am trying to solve the following differential equations: \begin{align} \dot{a}(\tau) &= \frac{a(\tau)}{\sqrt{3}}\left[ \rho_R(\tau) + \rho_V(\tau) \right]^{1/2}\\[0.25cm] \dot{\rho}_R(\tau) &+ 4\frac{\dot{a}(\tau)}{a(\tau)}\rho_R(\tau)=-\dot{\rho}_V(\tau) \end{align} with the initial conditions $$ a(\tau = 0) = 1,\qquad \rho_R(\tau = 0) = 1. $$ The function $\rho_V(\tau)$ is the source of all of the problems. It is of the form $\rho_V(r) = \exp(-I(\tau))$, where $$ I(\tau) = \frac{4\pi}{3} \int_{0}^{\tau}\Gamma(\tau')a^3(\tau')R^3(\tau', \tau)\ d\tau' $$ and $$ R(\tau', \tau) = \int_{\tau'}^\tau \frac{d\tau''}{a(\tau'')}, \qquad \Gamma(\tau) = \Gamma_0e^{\beta\tau} $$ Let $\Gamma_0 = \beta = 1$ for simplicity.
I am having trouble integrating this code using NDSolve since the integrals depend on $a(\tau)$. I would imagine that this is possible to evaluate since the integration only occurs until the currently evaluated time-step, but I am unsure how to do so. Below is my best attempt, but it doesn't work.
The key issue, is that I am unsure how to pass $a(\tau)$ into the function Ii so that I can perform the integration. Thanks for the help!
Γ[t_] := Γ0*Exp[β*t] ; Γ0 = 1; β = 1;
Rr[𝜏p_?NumericQ, 𝜏_?NumericQ, A_?NumericQ] := NIntegrate[1/A, {t, 𝜏p, 𝜏}];
Ii[𝜏_?NumericQ, A_] := (4*Pi)/3 NIntegrate[Γ[𝜏p]*A^3 Rr[𝜏p, 𝜏, A]^3, {𝜏p, 𝜏0, 𝜏}];
𝜌V[𝜏_ , A_?NumericQ] = Exp[-Ii[𝜏, A]];
𝜏0 = 0; 𝜏f = 50;
eqs = {
a'[𝜏] == a[𝜏]/Sqrt[3]*Sqrt[𝜌R[𝜏] + 𝜌V[𝜏, a[𝜏]]],
𝜌R'[𝜏] + 4*a'[𝜏]/a[𝜏] 𝜌R[𝜏] == -D[𝜌V[𝜏, a[𝜏]], 𝜏],
𝜌R[𝜏0] == 1, a[𝜏0] == 1
};
sols = NDSolve[eqs, {𝜏, 𝜏0, 𝜏f}];
Update 1
@Alex Trounev has seemed to produce many solutions, all of which agree with each other. I tried to develop another solution, which seems to get close but isn't quite the same as the three results below.
Much like, @Alex Trounev I defined
\begin{equation} R(\tau) = \int_0^\tau \frac{d \tau'}{a(\tau')}, \qquad \dot{R}(\tau) = \frac{1}{a(\tau)} . \end{equation}
Given $I(\tau)$ above, you can expand the cubic term and then define four variables
\begin{equation} v_i(\tau) = \int_{0}^{\tau}\Gamma(\tau')a^3(\tau) R^i(\tau)\ d\tau \end{equation}
each of which corresponds to the differential equation
\begin{equation} \dot{v}_i(\tau) = \Gamma(\tau')a^3(\tau) R^i(\tau), \qquad i = 0,1,2,3 \end{equation}
This allows us to write the function $I(\tau)$ as
\begin{equation} I(\tau) = \frac{4\pi}{3}\left\{ R^3(\tau)v_0(\tau) - 3R^2(\tau)v_1(\tau) + 3R(\tau)v_2(\tau) + v_3(\tau) \right\} . \end{equation}
The initial conditions for each $v_i(\tau = 0) = 0$. I tried integrating this system (see below). While I get a similar result, the height of the peak is significantly different. It also seems like there is an issue with Gam0 = 1. I'm not particularly sure what is causing the cusp at t ~ 2.
I am not sure where my error is here, so any input would be appreciated.
Gam[t_] := Gam0*Exp[beta*t] ; Gam0 = 1; beta = 1;
t0 = 0; tf = 4;
Ii[t_] := (4 Pi)/
3 (Rr[t]^3 v0[t] - 3* Rr[t]^2 v1[t] + 3* Rr[t]*v2[t] - v3[t]);
eqs = {
a'[t] == a[t]/Sqrt[3]*Sqrt[rhoR[t] + Exp[-Ii[t]]],
rhoR'[t] + 4*a'[t]/a[t] rhoR[t] == -D[Exp[-Ii[t]], t],
Rr'[t] == 1/a[t],
v0'[t] == Gam[t]*a[t]^3,
v1'[t] == Gam[t]*a[t]^3*Rr[t],
v2'[t] == Gam[t]*a[t]^3*Rr[t]^2,
v3'[t] == Gam[t]*a[t]^3*Rr[t]^3,
rhoR[t0] == 1,
a[t0] == 1,
Rr[t0] == 0,
v0[t0] == 0,
v1[t0] == 0,
v2[t0] == 0,
v3[t0] == 0
};
sols = NDSolve[eqs, {a, rhoR, Rr, v0, v1, v2, v3}, {t, t0, tf},
WorkingPrecision -> MachinePrecision];
Update 2
My issues was one rouge minus sign in front of v3 (corrected above). The code now reproduces @Alex Trounev's quite well. Problem solved!
