I've been working with Mathematica for symbolic calculations, primarily Integrations. The following is a function I am trying to integrate with the shown assumptions:

Integrate[(50 (0.000258748 + 0.117362 Vch^2) (23 + (84 Vch)/Sqrt[
   1 + 776.46 Vch^2]) (1 + (
   0.76354 (2 + 2329.38 Vch^2 + 1.20578*10^6 Vch^4))/(1 + 
     776.46 Vch^2)^(3/2)))/(1 + Vch^2), Vch, Assumptions -> {Vch != 0, Element[{Vch}, Reals]}]

However, I am getting various results for the output of the function depending on multiple criteria:

1)If I expand the above function before applying it to Integrate[], I get 1/0^2 infinity errors.

2) If I do not expand the function and leave it as shown above, the result are complex numbers

3) If I expand the above function and manually integrate each term (since the terms are added to each other) I will get a real result, which is what it should be.

I have to take multiple integrals similar to the one above. Taking each term after the expansions and Integrating each one is fine, but I think it's a roundabout solution.

Does anyone have any ideas as to what is happening and what steps I could take to get a real solution?

Thanks!