A positive hardware float (hfloat)  is "normally" in the range  10^(-308) .. 10^(308) [approx].
But  compromising precision, a subnormal representation (having the mantissa starting with zeros) allows smaller values up to about 5*10^(-324).
Your number 9.395*10^(-315) requires such a subnormal hfloat.
Note that evalhf  does conversions float -> hfloat -> float.
It seems that there is a bug in these conversions for subnormal hfloats. The bug seems to be very old.

restart;
interface(prettyprint=0):
Digits:=20:
evalhf(9.87654321*10^(-315));
-0.85702008631607784918e-314
evalhf( proc()9.87654321*10^(-315)end ());
0.
9.87654321*10^(-315); evalhf(%);
0.98765432100000000000e-314
-0.85702008631607784918e-314
k:=315: evalhf(9.87654321*10^(-k));
0.98765431971986665421e-314
for n from 305 to 330 do
  evalhf(9.87654321*10^(-n));
od;
0.987654321000016695e-304
0.987654321000016829e-305
0.987654321000016871e-306
0.987654321000017055e-307
0.98765432100001643399e-308
0.98765432099999568324e-309
0.98765432099994133602e-310
0.98765432099850854564e-311
0.98765432101431864631e-312
0.98765432095503076881e-313
0.98765431971986665421e-314
0.98765430588602857065e-315
0.98765456774082086651e-316
0.98765303613731875865e-317
0.98764216669311025122e-318
0.98763722603665183876e-319
0.98566096345328685259e-320
0.97824997876566815445e-321
0.98813129168249308530e-322
0.
0.
0.
0.
0.
0.
0.

At each step except the first, you take tau12, tau21 from the previous step.

Numeric only (of course).

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(2)=1, D(y)(2)=1/2], numeric):
plots:-odeplot(ds, r=1e-6..100, color=red);

Nothing unusual. The integrator is very complex, so, equivalent functions could be integrated using different methods.

But it's easy to have the same unevaluated result using:
int(simplify(integrand2), x);  ##  :-)

In the absence of a better solution (the generated LaTeX code is not easy to manipulate), I'd replace in the worksheet e.g.

sys_set := convert(sys_vec, set);

with

sys_set := convert(sys_vec, set):
for eq in sys_set do print(eq) od; # display elements

You don't need Maple for this: z / |z|^2 = z is holomorphic, so the answer is 0.

The limit of a function F : A --> B, at a point u is defined when u is an accumulation point (= limit point = cluster point) of A.
In our case the standard choises for A, B are B = R and A = R \ {0} \ {1/(n*Pi) : n in Z\{0}}.
0 is an accumulation point for A, so the limit problem makes sense. It exists and equals 1 (using the standard limit : sin(t)/t, for t --> 0).

[A more techical note: the topology in A is by default the subspace topology induced by the natural topology in R; your neighbourhood refers to this topology]. 

Your ode is

ode := diff(y(x), x) = 15*exp(4*y(x)) - 5*y(x)

There is no explicit solution for this ODE; this is why the Qualitative theory of ODEs exists.
You can study the solutions taking y0 = y(0) as a parameter. You will see that the solution always blows up in a finite interval.

Actually in our case it's easy to see that f(y)>0 (=  rhs(ode)) and

int(1/f(y), y=y0 .. y(x)) = x - 0

so, the existence domain for x>0 is bounded, because LHS is obviously bounded.

P.S. It is a good idea to post the problem, not just in the worksheet; it is short, and the worksheet should be only optional. You will increase the chances to get answers.

Your function has an essential singularity at z=0. Maple (and probably other CAS too) cannot compute the series for such singularities. (The series has the form Sum(a(n)*z^n, n = -infinity .. infinity)).

But you can use Maple to obtain it by multiplying two series (try it first by hand).

This occurs frequently when the symbolic solution is expressed via mathematical functions with branches.

The numerical solution is continuous (in general), and it may coorespond to different branches (i.e. it "jumps" from one branch to another).

``||i   are valid names, so the assignments work.

But note that they are global variables. So, if you have

local a, b__c, t;

then the assignments will produce new global variables, e.g. :-a  will be 1, but a remains free (unassigned).

value(IntegrationTools:-Change(z, s=x-t, t)) assuming x>0;

But do you really need it?

You simply forgot to load the package. Start with:

with(Student[Calculus1]):

Selecting the homogeneous solutions is not complicated:

with(PDETools):
PDE := 2*((x1 - x2)^2 + (y1 - y2)^2)*(u1*diff(f(x1, y1, x2, y2, u1, v1, u2, v2), x1) + v1*diff(f(x1, y1, x2, y2, u1, v1, u2, v2), y1) + u2*diff(f(x1, y1, x2, y2, u1, v1, u2, v2), x2) + v2*diff(f(x1, y1, x2, y2, u1, v1, u2, v2), y2)) - ((u1 - v1)^2 + (u2 - v2)^2)*((x1 - x2)*diff(f(x1, y1, x2, y2, u1, v1, u2, v2), u1) + (y1 - y2)*diff(f(x1, y1, x2, y2, u1, v1, u2, v2), v1) - (x1 - x2)*diff(f(x1, y1, x2, y2, u1, v1, u2, v2), u2) - (y1 - y2)*diff(f(x1, y1, x2, y2, u1, v1, u2, v2), v2)) = 0:
n:=2:
P:=rhs(PolynomialSolutions(PDE, degree = n)[]):
C:=[indets(P,suffixed(_C, integer))[]]:
T:=[coeffs(collect(P,C,'distributed'),C)]:
select(u -> degree(u)=n, T); 

[v1^2 + 2*v1*v2 + v2^2, (v1 + v2)*u1 + (v1 + v2)*u2, u1^2 + 2*u1*u2 + u2^2, u1*y2 + u2*y1 - v1*x2 - v2*x1, u1*y1 + u2*y2 - v1*x1 - v2*x2]

But it seems that you want to eliminate also other terms. This will need extra work (probably using Groebner bases).

There are essentially two methods.

1. Use DirectSearch.
GlobalOptima(f(x,y), {x=3..4, y=0..1}, maximize);
The approx max is (almost) 0.

2. Use calculus.
Solve the system [diff(f(x,y),x), diff(f(x,y),y)].
Maple finds a unique solution in the domain: x=4, y=0, ( f(4,0)=0 ).
It remains to compute the max on the boundary: f(x,0), f(x,1), f(3,y), f(4,y)
(it's easy).