@Carl Love 

My point is that evalf/ln1 should use a special algorithm for ln(1+x) (e.g. via series) being so much faster. BTW, Digits=64 in ln1 is enough for OP's example.

@Carl Love 

Note that this way in OP's example the procedure evalf/ln1 uses Digits=91 (!). So, what it does is:

evalf[32](evalf[91](ln(1+exp(-64))));

All these symbols are defined here: https://en.wikipedia.org/wiki/List_of_mathematical_symbols

It's a matter of preferences and definitions/notations, rather than ambiguity. Even N could represent {0,1,...} or {1,2,...}, but the exact definition is clear from the context, or reminded in the case of a serious textbook or article.

@Carl Love 

It is unfortunate that they are not documented.

restart;
f := x -> sin(2*x)+3;
g := x -> sin(3*x)+6;
a:=0; b:=5;
plots:-animate(plot3d,
[ [ [a+t*(b-a), f(a+t*(b-a))*cos(u), f(a+t*(b-a))*sin(u)], 
    [a+t*(b-a), g(a+t*(b-a))*cos(u), g(a+t*(b-a))*sin(u)],
    [a, (f(a)+t*(g(a)-f(a)))*cos(u), (f(a)+t*(g(a)-f(a)))*sin(u)],
    [b, (f(b)+t*(g(b)-f(b)))*cos(u), (f(b)+t*(g(b)-f(b)))*sin(u)]
  ], 
t=0..1, u=0..U, transparency=0.3, grid=[40, max(2,floor(7*U))]
],  U = 0..2*Pi, axes=none,  paraminfo=false);

@Carl Love 

Please note that I have mentioned the sfloat->hfloat->sfloat conversion from the beginning. My impression was that the OP obtained the result in a evalhf environment and it was distorted. But now I see that he used a higher precision and wants to save the result as hfloat. For this of course it's enough to save it as float[8] without any evalhf or option hfloat.

@AOdrzywolek 

You are wrong here. It is proved mathematically that a long double can be exactly retrived (each of its 64 bits) using 18 decimal digits (it is irrelevant how it was obtained). [ The fact that to obtain n correct digits with a math function you may need Digits>n is quite another thing: after the computation you save evalf[18](result) ].

@AOdrzywolek 

You can use Digits=18. Then it seems the bug does not occur. But for portability, Carl's advise for ASCII is a better solution (even if the file could be up to 3 times bigger).

@Carl Love 

Sorry, but I cannot agree. I still think that it's a bug; probably even the the 20 digits display for subnormal hfloats is not intended.

Why should we decrease Digits when using evalhf? For "regular" hfloats, everything is OK. 

@Carl Love 

Note that evalhf displays 18 digits. That's because this way the exact hfloat representation can be retrived.

Digits:=10:
evalhf(123456.7890123456789);

        123456.789012345675

I don't see any contradiction using evalhf when Digits is 30. It could be e.g. a partial result, or a first approximation.

@Axel Vogt   Sorry, I don't know.

@zjervoj97 

Probably you must define tau12, tau21 before  lng1 := ...
I's just a guess, you should see what quantities change at each step. Only you know that.

@SIMBA 

y = 0 obviously satisfies your conditions. Probably you want something else, but then try to impose conditions for uniqueness (if you are sure about existence; I am not).

Here is an Explore procedure for a visual "shooting":

a:=1e-1:
ode:=((D@@2)(y))(r)+(D(y))(r)/r-sin(y(r))*cos(y(r))/r^2+sin(y(r))^2/r-sin(y(r))-sin(2*y(r)) = 0:
ds:=dsolve([ode, y(a)=A, D(y)(a)=B], numeric, parameters=[A,B]):
P:=proc(A1,B1) 
   ds(parameters=[A=A1,B=B1]); 
   plots:-odeplot(ds, [r, y(r)], r = a .. 100, view=-10..10);
end proc:
Explore( P(A1,B1),parameters=[A1=-1.0 .. 1.0, B1=-5.0 .. 5.0] );

@SIMBA 

Ordinary differential equations (ODEs)  are studied in most scientific faculties in universities. There are many textbooks. I will not recommend any of them because I do not know your background.
Maple has two main approaches to solve ODEs: symbolically  (exact) and numerically. Unfortunately exact solutions exist only for very few equations/systems.
The main Maple procedure for solving ODEs is dsolve. Read as much as you can in the help system. Be aware that the quantity of information is high, but there are also many examples.

About your last question. ds is a numeric procedure returned by dsolve. For example,

ds(3.0);
returns
[r = 3., y(r) = 1.91101863645892, diff(y(r), r) = 1.03474512991423]

i.e. it gives the numerical value of the dependent variable and its derivative at r = 3.0.

@Carl Love 

Actually, when I said "the result is correct" I meant "I cannot find a counterexample". You changed everywhere ln(exp(u)) with u; this is not always true.

Edit. I found one:

f:= z-> ln(exp(10/z - 10/sqrt(z^2))):

@Carl Love 

OK, now the result is correct but it's very far from your usual style. Instead of a direct and simple formula, you prefer Change and then manipulate the result. You could do this as well:

r:= t-> exp(I*t):
call_anything_such_as_Change():
subs(%=Int(f(r(t))*r(t), t=-%Pi.. %Pi)*I, %);