A quick question, does first order ODE system always have numerical solutions?

I mean like, just say a 2 ODE system, there are parameters A=100, B=10 and it has numerical solution y(1) and y(2) for the given initial time and final time.
Now is it possible that if A = -100, B=10 then it does not have correct numerical solution?

When such issue happens,
basically the Fortran ODE solver we are using will give some error status, right?
Or the ODE solver just crash or something?

That really depends. In that respect a numerical solver is quite different than the mathematics. You would hope that a solver is capable of determining that. Here is a very simple ODE:

dy/dx = y**2 + 1, y(0) = 0

The matematical solution is y = tan(x) (if I have not made a mistake), so there is a singularity at x = pi/2. This is no particular problem for the mathematics, but the numerical solution can be dramatically wrong around the singularity and it may not be immediately clear that that is the case!

I agree with Arjen that there is a distinction between the numerics and the mathematics. Disclaimer, I havenâ€™t used the Fortran ODE solver:

Mathematics part: you could cook up an ODE that does not have a solution or whose solution changes dramatically when you change the parameters of the problem. So in theory you need well-posedness of the problem.

Numerics part: the equations of your scheme may be ill-conditioned, meaning they are sound mathematically, but the finite precision of computers poses a problem, like in the case of Arjenâ€™s example. Here you can expect most solvers to warn you.

For PDEs, you can cook up problems that are well-posed, the final matrix can be well-conditioned, but it is just the wrong numerical scheme for the problem: youâ€™ll obtain a wrong solution without any warnings. I am less familiar with ODEs.

I too agree with Arjen, and would add some jargon for you to watch out for: â€śstiffâ€ť equations. These are where some of the derivatives have wildly different absolute values to others. So for some equations the step length is far too small, and rounding error is obscuring real processes in the maths, and the same step length is far to big for other equations in the system and the error in the numerical approximation has become unacceptably too large. In these cases it is worth looking at Gear-Hindmarsh methods and their more modern successors.

Uhm, to be more specific, I am using @prince_mahajan 's FLINT ODE solver,

It is a fast ODE solver with adaptive steps.
But in some cases, I found that, even if it reaches the max steps (using its default DOP853 method), it still cannot complete some ODE integrals. Increasing the max steps does not help.

So I am thinking how to finish the numerical integrals anyway. Perhaps I can try different methods other than DOP853 etc.

Or, I am just wondering, is it just that, some ODEs just do not have numerical solutions within the given time range and initial condition? If so no matter what ODE solver or method I use, perhaps it just do not have a solution.

I have not used FLINT, but it appears to be Runge Kutta and explicit. So you may want to look at even quite simple implicit methods which can be unconditionally stable.(But potentially expensive)
The other route to go are the predictor-corrector methods, the most famous of which are the Adams-Bashforth schemes, if the RK stiffness detector is issuing a warning then go to Gear-Hindmarsh.