# Is it possible to solve the following system of integrodifferential equations in Mathematica?

Consider a system of equations
$$begin{cases}F(E,t) = frac{partial f(E,t)}{partial t}-H(T_{gamma})Efrac{partial f(E,t)}{partial E} -I(f(E,t),E,T_{gamma})=0, \ G(t) = frac{partial T_{gamma}}{partial t} + H(T_{gamma})T_{gamma}+int limits_{0}^{infty} dE I(f(E,t),E,T_{gamma})=0,end{cases}$$
where
$$I(f(E,t),E,T_{gamma}) equivint limits_{Gamma(E)} dE’ G_{1}(E’,E,f(E,t),f(E’,t),T_{gamma})$$
is some integral and $$H$$ is a known function and $$Gamma(E)$$ is the domain of integration.

Could you please tell me whether it is possible to solve this system in Mathematica? As far as I know, NDSolve does not work with integrodifferential systems of equations, but maybe some scheme exists?

These equations in Mathematica code, `Equationfn(En, t)`, `EquationTg(t)` are given below (there is some piece of code that defines functions):

``````ME = 0.5;
etaBPlanck = 6.1*10^-10;
ng(T_) = 2*Zeta(3)/Pi^2 T^3;
rhoeTemp(T_) :=
4/(2*Pi^2)
NIntegrate((Ee^2 Sqrt(Ee^2 - ME^2))/(
Exp(Ee/T) + 1), {Ee, ME, Infinity})
peTemp(T_) :=
4/(6*Pi^2)
NIntegrate((Ee^2 - ME^2)^(3/2)/(Exp(Ee/T) + 1), {Ee, ME, Infinity})
rhoe(T_) =
Quiet(10^Interpolation(
Join(Table({10^T, Log10(etaBPlanck*ng(10^T)*ME)}, {T,
Log10(10^-8), Log10(0.0049), 0.1}),
Table({10^T,
Log10(etaBPlanck*ng(10^T)*ME + rhoeTemp(10^T))}, {T,
Log10(0.005), Log10(19.99), 0.01}),
Table({T, Log10(7/8*4*Pi^2/30*T^4)}, {T, 20, 300, 0.5})),
InterpolationOrder -> 1)(T))
pe(T_) = Quiet(
10^Interpolation(
Join(Table({10^T, -90}, {T, Log10(10^-8), Log10(0.0049), 0.1}),
Table({10^T, Log10(10^-99 + peTemp(10^T))}, {T, Log10(0.005),
Log10(19.99), 0.01}),
Table({T, Log10(1/3 7/8*4*Pi^2/30*T^4)}, {T, 20, 300, 0.5})),
InterpolationOrder -> 1)(T))
gStare(T_) := rhoe(T)/(Pi^2/30 T^4)
DrhoeDT(T_) = D(rhoe(T), T);
(*Neutrino and photon*)
rhon(T_) = 7/8*6*Pi^2/30*T^4;
rhona(T_) = 7/8*2*Pi^2/30 T^4;
rhog(T_) = 2*Pi^2/30 T^4;
DrhogDT(T_) = D(rhog(T), T);

GFvalue = 1.16*10^-11;
glValue = 1/2 + 0.22;
grValue = -(1/2) + 0.22;
hbar = 6.582119*10^-25;
sToGeVminusOne = 1/hbar;
sToMeVminusOne = sToGeVminusOne*10^-3;
CollisionIntegralnnToee(En_,
t_) := -fn(En, t)*
NIntegrate(
sToMeVminusOne (
8 En E2^3 GFvalue^2 (glValue^2 + grValue^2) fn(E2, t))/(
9 (Pi)^3), {E2, 0, Infinity}) + (
16 E^(-(En/Tg(t))) En *sToMeVminusOne*
GFvalue^2 (glValue^2 + grValue^2) Tg(t)^4)/(3 (Pi)^3)
CollisionIntegralneTone(En_,
t_) := -((
64 En*sToMeVminusOne*
GFvalue^2 (glValue^2 + grValue^2) Tg(t)^4 fn(En, t))/(
3 (Pi)^3)) +
NIntegrate(
sToMeVminusOne/(En^2 (Pi)^3) 8 E^(-(En/Tg(t)))
GFvalue^2 (glValue^2 + grValue^2) Tg(
t)^2 (-En^2 (E3^2 + 2 E3 Tg(t) -
2 (-1 + E^(E3/Tg(t))) Tg(t)^2) -
2 En Tg(t) (E3^2 + 2 E3 Tg(t) -
2 (-1 + E^(E3/Tg(t))) Tg(t)^2) +
2 Tg(t)^2 ((-1 + E^(E3/Tg(t))) E3^2 -
2 (1 + E^(E3/Tg(t))) E3 Tg(t) +
4 (-1 + E^(E3/Tg(t))) Tg(t)^2)) fn(E3, t), {E3, 0, En}) +
NIntegrate(
sToMeVminusOne/(En^2 (Pi)^3) 8 E^(-(En/Tg(t)))
GFvalue^2 (glValue^2 + grValue^2) Tg(
t)^2 (2 (-1 + E^(En/Tg(t))) Tg(
t)^2 (E3^2 + 2 E3 Tg(t) + 4 Tg(t)^2) -
En^2 (E3^2 + 2 E3 Tg(t) - 2 (-1 + E^(En/Tg(t))) Tg(t)^2) -
2 En Tg(t) (E3^2 + 2 E3 Tg(t) +
2 (1 + E^(En/Tg(t))) Tg(t)^2)) fn(E3, t), {E3, En, Infinity})
H(t_) = sToMeVminusOne (1.66 Sqrt(10.75) Tg(t)^2)/(1.2*10^22);
Equationfn(En_, t_) := ((D(fn(En, tt), tt)) /. {tt -> t}) -
H(Tg(t)) En*((D(fn(en, t), en)) /. {en -> En}) -
CollisionIntegralnnToee(En, t) - CollisionIntegralneTone(En, t)==0
EquationTg(t_) := ((D(Tg(tt), tt)) /. {tt -> t}) +
1/(DrhoeDT(Tg(t)) +
DrhogDT(Tg(t))) (4*H(Tg(t))*rhog(Tg(t)) +
3*H(Tg(t))*(rhoe(Tg(t)) + pe(Tg(t))) +
NIntegrate(
CollisionIntegralnnToee(En, t) En^3/(2*Pi^2), {En, 0,
Infinity}) +
NIntegrate(
CollisionIntegralneTone(En, t) En^3/(2*Pi^2), {En, 0,
Infinity}))==0
``````