Just run the model with solvers from 3 languages: Julia, C++, and Fortran.
using OrdinaryDiffEq, Sundials, ODEInterfaceDiffEq
sol1 = solve(prob1,Rodas5(),abstol=1e-12,reltol=1e-12)
sol2 = solve(prob1,CVODE_BDF(),abstol=1e-12,reltol=1e-12)
sol3 = solve(prob1,radau(),abstol=1e-12,reltol=1e-12)
@show sol1.t[end], sol2.t[end], sol3.t[end]
@show sol1[end], sol2[end], sol3[end]
(sol1.t[end], sol2.t[end], sol3.t[end]) = (14.57381167838683, 14.57381163596623, 14.573811679098052)
(sol1[end], sol2[end], sol3[end]) = ([213178.70363115682, 4.8383636378845696e17], [1.868111482362633e36, 2.930079193719137e117], [349223.9235893567, 2.3765960233233695e18])
You see that every robust solver from each language all blows up to infinity at almost precisely the same point in time. That is very strong evidence that the true solution of the model you had written down goes to infinity at that time.
It seems like the automatic stiffness detection falls into a loop (which we should fix), but at least your issue is that the differential equation you wrote down is divergent and you should probably double check for sign issues.