@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, %);

@Carl Love 

The simplify/symbolic may affect f because it could contain ln(exp(...)).
[A tested counterexample is f:= z-> ln(exp(1/z - 1/sqrt(z^2))) ]
You try do obtain: Int(f(exp(t*I))*exp(t*I), t = -%Pi .. %Pi)*I
Why not use it directly?

@Carl Love 

I don't find a good idea to use Change in this case. Trying to repair with simplify/symbolic has no chance.
For example, f := z -> 1/z - 1/sqrt(z^2)

@Carl Love 

Your result will be 0 for any f.

@Kitonum 

The other solution is also traditional: integration by parts.

@Carl Love 

Probably simplify just gives 0 for Sum(0,...). But it's nice to see the huge amount of work Maple does to simplify (using the assumptions on n): use 
infolevel[all]:=5

@mmcdara 

If the variables are already assigned, you will have to change a bit the code:

restart;
vars   := [a, b__c, t[1]];
a:=1234;
assign(eval(vars,1) =~ [$1..3]);

@acer 

A useful info!

It works in Windows 64 bit if  ':-LIB'="libmkl_rt.so"  is replaced by ':-LIB'="mkl_rt"