@acer 

Of course, if sigma is very small, increasing Digits will be mandatory.

@sand15athome 

You have two solutions: switch to software floats or simplify the function using Heaviside (in the latter case evalhf has no problems).
You can inspect `f*g`(x). You will see that with my sugestion it does not contain Dirac and is much simpler (and probably also faster).

So, you have a continuous function Phi(t). It seems to be (almost) periodic.
What do you mean by a probability density associated to it?

@sideshow 

The solution is not unique. E.g.

Sol := proc (T) options operator, arrow; fsolve({FM(q, Jo, T) = 1, FQ(q, Jo, T) = q}, {Jo=0.5 .. 1.5, q=0.5 .. 1}) end proc:Sol(0.3);
              {Jo = 1.124786970, q = 0.7332828274}

@leafgreen 

Your function f(x) = 1/ln(x + 1/x)  does not have an elementary antiderivative (= indefinite integral). This can be proved using https://en.wikipedia.org/wiki/Risch_algorithm.
Maple applies it for this kind of functions. Unfortunately such algorithm does not exist for definite integrals.

Most such integrals cannot be expressed with elementary (or even special) functions.
Why do you think that this one must be an exception?

@sideshow 

fsolve gives  {Jo = -1.082524999, T = -.4330099990}.

What does not work well?

@sideshow 

I have considered the system

FM(q, Jo, T) = 1,   FQ(q, Jo, T) = 2.

That's because the original
FM(q, Jo, T) = 0,   FQ(q, Jo, T) = 0,
obviously has no solution. And the Mathematica code contains something similar.
You must set the system you want to solve, I cannot guess it.

@Jonas123456 

Just curious. What is the source of your expression?

@acer 

I did not do it by hand, but Maxima agrees.

Edit. Maybe  99.9999 because Derive agrees too (yes, I still use it sometimes!)

This reminds me of:

@williamdiao 

The only problem is that Quo is called with 2 or 3 arguments. So, e.g. use:

C[Quo] := proc(a,b,r) if nargs=3 then quo(a,b,z,r) else quo(a,b,z) fi end:

@Math111 

Then what is (should be) F''(0) ?
[the fact that a method provides a solution does not imply that it actually exists].

Edit.

RK-4 for a BC problem? Probably some shooting method.

The ODE has the following explicit form:

The denominator has big chances to be 0 for r in [0,1].

Note that:
|F'(r)| must be < 0.25 ==>  |F(r)| < 1 + 0.25 r

The usual way of presenting this would be:

Write the projection as a procedure. E.g. the projection onto the unit sphere is
P:=(x,y,z) -> (x/sqrt(x^2+y^2+z^2), y/sqrt(x^2+y^2+z^2), z/sqrt(x^2+y^2+z^2));

Define the curve. E.g.
C:=t -> (cos(t),sin(t),t);

Plot the projection:
plots:-spacecurve( [P(C(t))],t=0..8*Pi);
or
plots:-animate(plots:-spacecurve, [[P(C(t))],t=0..a],a=0..8*Pi );

Then your posts will have a larger audience.

@Simwar 

For rather small intervals you can check all the possible cases.

Once you have a plausible formula, you could try to compute the sum using Maple (with assume facility, maybe fixing some parameters). Maple is "competent" in computing this kind of sums.

But if you have derived correctly the formula, a comprehensive random test should be enough for your peace of mind.