How to define an appropriate upper bound integral function

Question

Recently, I am doing a course project, and the professor forces us to use Mathematica; however, I am new to Mathematica, and here are my problems.

Setup

I am required to, eventually, numerically solve a group of differential equations which look like the form: $$\frac{dX(t)}{dt}=P(X(t))+F(t) $$ Here, $X(t)$ is the coloumn vector $(x_1(t),x_2(t),x_3(t))^T$, $P(X(t))$ is a polynomial of $X(t)$ which contains terms like $x_i(t)x_j(t)$.

$F(t)$ is the key term for this project. $F(t)$ is a column vector takes form $(F_1[x_1(t)]\cdot x_1^2(t),\,F_2[x_2(t)]\cdot x_2^2(t),\,F_3[x_3(t)]\cdot x_3^2(t))$, where $F_i[x_i(t)]$ is a upper bound function of a complicated function.

Here, I will just present $F_1[x_1(t)]$, since $F_2$ and $F_3$ are similar to $F_1$. I think I can deal with $F_2$ and $F_3$ once I know the way to deal with $F_1$.

$$F_1[x_1(t)]:=\int_0^{\frac{x_1(t)}{3}}\frac{(1+y)(2y-1)^2(2y-3)^2\sin^3(\pi y)[\Gamma(y-1)]^2\Gamma(-2y)}{(y-2)\pi^3}dy $$ Here $\Gamma$ is the usual gamma function. So far, we are given that $F_1(x)$ has poles at $x=\frac{15}{2}+3n$. Similarly to $F_2$ and $F_3$. Here, we only care about the value before the first pole.

Main Problem

Now, here is my question. In order to plug $F_1[x_1(t)]$ in to the differential equation, I tried to numerically compute $F_1(x)$ and define $F_1(x)$ as the following code:

I1[x_] := (1 + x)*(2*x - 1)^2*(2*x - 3)^2*Sin[(\[Pi]*x)]^3*
   Gamma[(x - 1)]^2 Gamma[(-2)*x]/((x - 2)*\[Pi]^3)

F1[x_] := NIntegrate[I1[t]*Boole[t < 2.5 - Exp[-NF2]], {t, 0, x/3}]

Since we only care $0$ to $15/2$, I set a Boole Function take the value from $0$ to $2.5-\text{Exp(-NF2)}$ $[\text{Note: }15/2\div 3]$, where NF2 is some large number which is given by professor.

After I define the $F_1(x)$, I plug the whole $F_1[x_1(t)]$ into the differential equation, then I get a lot of warning. Can someone help me to fix the issue? Thanks.

Attached is the graph of $F_1(x)$ and the error for directly plug-in procedure. $\alpha(t)$ is $x(t)$ here.

Code

Here is the code I am working on.

NF3 = 40; NF2 = 24; Y2 = 
 1/2; Y3 = 0; DR22 = 2; DR32 = 1; DR23 = 1; DR33 = 3; \[Alpha]y = 0;
b1 = 41/3 + 8/3*Y2^2*NF2*DR22*DR32;
b2 = 19/3 - 4*NF2*DR32/3;
b3 = 14 - 4*NF3*DR23/3;
c1 = 199/9 + 8/3*Y2^4*NF2*DR22*DR32;
c2 = 35/3 + 49*NF2*DR32/3;
c3 = -52 + 76*NF3*DR23/3;
d1 = 88/3 + 32/3*Y2^2*NF2*DR22*DR32;
d2 = 24 + 16/3*NF2*DR32;
d3 = 9 + 3*NF3*DR23;
e1 = 9 + 6*Y2^2*NF2*DR22*DR32;
e2 = 3 + 4*Y2^2*NF2*DR32;
e3 = 11/3 + 4*NF3*Y3^2*DR23;

NC2 = 2;
NC3 = 3;
I1[x_] := (1 + x)*(2*x - 1)^2*(2*x - 3)^2*Sin[(\[Pi]*x)]^3*
   Gamma[(x - 1)]^2 Gamma[(-2)*x]/((x - 2)*\[Pi]^3);
F1[x_] := NIntegrate[I1[t]*Boole[t < 2.5 - Exp[-NF2]], {t, 0, x/3}];
I22[x_] := (NC2^2 - 1)/NC2 + (20 - 43 x + 32 x^2 - 14 x^3 + 4 x^4)*
    NC2/(2*(2 x - 1)*(2 x - 3)*(1 - x^2));
I23[x_] := (NC3^2 - 1)/NC3 + (20 - 43 x + 32 x^2 - 14 x^3 + 4 x^4)*
    NC3/(2*(2 x - 1)*(2 x - 3)*(1 - x^2));
H12[x_] := -11/2*NC2 + 
   NIntegrate[I1[t]*I22[t]*Boole[t < 1 - Exp[-NF3]], {t, 0, x/3}];
H13[x_] := -11/2*NC3 + 
   NIntegrate[I1[t]*I22[t]*Boole[t < 1 - Exp[-NF3]], {t, 0, x/3}];
Plot[F1[t], {t, 0, 8}]

A1 = 4*\[Alpha]1[t]*NF2*Y2^2*DR22*DR32;
A2 = 2*\[Alpha]2[t]*NF2*DR32;
A3 = 2*\[Alpha]3[t]*NF3*DR23;

odeg1 = 2*A1*\[Alpha]1[t]*F1[A1]/(3*NF2) + 
    D[\[Alpha]1[t], 
     t] == (b1 + c1*\[Alpha]1[t] + d1*\[Alpha]3[t] + 
      e1*\[Alpha]2[t] - 17*\[Alpha]y/3)*\[Alpha]1[t]^2;
odeg2 = 2*A2*\[Alpha]2[t]*H12[A2]/(3*NF2) + 
    D[\[Alpha]2[t], 
     t] == (-b2 + c2*\[Alpha]2[t] + d2*\[Alpha]3[t] + 
      e2*\[Alpha]1[t] - 3*\[Alpha]y)*\[Alpha]2[t]^2;
odeg3 = 2*A3*\[Alpha]3[t]*H13[A3]/(3*NF3) + 
    D[\[Alpha]3[t], 
     t] == (-b3 + c3*\[Alpha]3[t] + d3*\[Alpha]2[t] + 
      e3*\[Alpha]1[t] - 4*\[Alpha]y)*\[Alpha]3[t]^2;
Ncorrect = 
  NDSolve[{odeg1, odeg2, 
    odeg3, \[Alpha]1[3] == 0.00084, \[Alpha]2[3] == 
     0.00256, \[Alpha]3[3] == 
     0.0061}, {\[Alpha]1, \[Alpha]2, \[Alpha]3}, {t, 0, 20}];
Plot[Evaluate[{\[Alpha]1[t], \[Alpha]2[t], \[Alpha]3[t]} /. 
   Ncorrect], {t, 0, 20}, PlotRange -> 1, AxesLabel -> {t, Strength}, 
 ImageSize -> Large]

Answer 1

The use of x_?NumericQ which means the function evaluates only for numerical values should fix your problem. However, to compute F1, instead of computing the (indefinite) numerical integral one can turn it into a first-order DE:

$f(x)=\int_0^x g(t) \, dt \implies f'(x)=g(x)$ with $f(0)=0$

and solve this using NDSolve instead, which is much faster as it computes the values of F1 over the entire range in one go. Also, the solution is returned as an InterpolatingFunction so the issue that you faced disappears. But, because of the singularities of $g(x)$ at $x=0$ and $x=5/2$, one has to deal with these separately. To do this, perform series expansion about $x=5/2$ and subtract of the resulting pole there, (which leads to a logarithmic singularity).