you would like to r

• rewrite a differential equation, let's say

• a*diff(f(x), x)+b*f(x)=c

in a symbolic form 

• (a*D + b*I)(f) = c

where I is the identity operator and D the differential one.

• Next to define the formal inverse L of (a*D + b*I)

  L = (a*D + b*I)-1
  

• And use L to simplify another ODE which contains f(x) by replacing each of its occurrences by L(c)

This is possible only in very specific situations (see at the end of the attached file)
This file unfolds the previous "program" step by step and shows its limits when placed in a general framework.

Operator_manipulations.mw

Thus, the answer to your question "can it be done by factoring differential operators and NOT directly using dsolve?" is unfortunatele NO (or it is far beyond my capabilities)

Identifying a list of integers to a random variable may be  necessary or not (for instance tandomly chose the rows or columns of a matrix, as it happens in some numerical algorithms).
As your list is not a list of contiguous or evenly space elements, maybe the statistical approach can be useful?

Here are a few alternatives based on an identification to a random variable (the last alternative doesn't).
So make your choice.
 

Download Alternatives.mw

One thing: you can set a=1 without any loss of generality (same if you set d=1).

With this in mind it becomes extremely easy to construct by hand (almost by hand) particular solutions.
In the attached file one respects the constraints that the coefficients are in [1, 10], another don't.

galimatias.mw

Of course it's very far from Kitonum's exhaustive search

Look here view.aspx (paragraph "The External C Compiler").

Contrary to Windows platforms
"On all UNIX and Mac OS X systems, the GNU C compiler is not distributed with Maple, and must be installed separately, if it is not already present on the system " 

This could be the cause of your problem

One way which introduces only a slight modification to your code is to replace the line

eq1 := Qbar = InvT . Q . R . T . InvR

by the red line below

assign([entries(Qbar, nolist)] =~ [entries(InvT . Q . R . T . InvR, nolist)]);

plot(Qb11, p = 0 .. 9): # fails for Qb11 depends on 3 variables:

Warning, expecting only range variable p in expression (.2339357430e12*cos(p)^2+.2319277108e11*s^2)*cos(p)^2+(.2319277108e11*cos(p)^2+s^2*Q22)*s^2+.2868e11*s^2*cos(p)^2 to be plotted but found names [Q22, s]

# check this
indets(Qb11, name);

                          {Q22, p, s}

# do this, for instance, to evaluate Qb11 for siome values of Q22 and s
plot(eval(Qb11, [Q22=1, s=1]), p=0..9)

with recent Maple versions

numer(a) %/ denom(a);

For older versions (2015 for instance)
 

with(InertForm):
n := convert(numer(a), string):
d := convert(denom(a), string):
Parse( cat("(", n, ")/", d) ):
Display(%);

Hints for Problem 7
An equivalence relation E has 3 properties R, S and T (which are?).

• one of them (R) is obvious
• one (S) comes from the commutativity of the addition
• the last (T) involes 3 couples C, C' and C"
  • write C E C' and C' E C""
  • try to combine the two equalities to obtain C E C" (needs a simplification and uses properties of the addition)
     
• In fact, let F a relation which has the same properties of the addition, then  (a, b) E (c, d) defined by  (a+d) F (b+c) is an equivalence relation.

For the three other problems look at some hints here
pbs.pdf

It's clear that D is at the intersection of the bisector of angle A and the line which passes through B and C.
Thus its coordinates and a parametric representation of AD

coords__D := ([-6, 1,-1] +~ [2,-3, 7]) /~ 2;
line__AD  := [x, y, z] =~ [-2, 3, 5]*~s +~ coords__D*~(1-s);

But as I'm fond of the geom3d package, here is a step by step wauy yo use it (uselesss for this problem but it could help you for more complex ones). 

restart:
with(geom3d):

point(A, -2, 3, 5):
point(B, -6, 1,-1):
point(C,  2,-3, 7):

line(AB, [A, B]):
line(AC, [A, C]):

# parametric equations of AB and AC

eq__AB := [x, y, z] =~ Equation(AB, s):
eq__AC := [x, y, z] =~ Equation(AC, s):

# AD is the bisector of angle BAC:
# parametric equation of AD

eq__AD := (eq__AB +~ eq__AC) /~ 2;
               [x = -2, y = 3 - 4 s, z = 5 - 2 s]

# implicit equation of AD 
map(u -> if depends(u, s) then solve(u, s) else rhs(u) end if, eq__AD):
eq__AD__imp := add(%);
                         5   1     1  
                         - - - y - - z
                         4   4     2  

# coordinates of D
local D:
line(BC, [B, C]):
line(AD, rhs~(eq__AD), s):
intersection(D, BC, AD):
detail(D)
          Geom3dDetail(["name of the object", D], 
            ["form of the object", point3d], 
            ["coordinates of the point", [-2, -1, 3]])

# visualize
plots:-display(
  PLOT3D(
    POINTS([coordinates(A)], SYMBOL(_SOLIDCIRCLE, 30), COLOR(RGB, 1, 0, 0)),
    POINTS([coordinates(B)], SYMBOL(_SOLIDCIRCLE, 30), COLOR(RGB, 1, 0, 0)),
    POINTS([coordinates(C)], SYMBOL(_SOLIDCIRCLE, 30), COLOR(RGB, 1, 0, 0)),
    POINTS([coordinates(D)], SYMBOL(_SOLIDCIRCLE, 30), COLOR(RGB, 0, 0, 1)),
    TEXT(1.1*~coordinates(A), "A", FONT(Helvetica, bold, 14)),
    TEXT(1.1*~coordinates(B), "B", FONT(Helvetica, bold, 14)),
    TEXT(1.1*~coordinates(C), "C", FONT(Helvetica, bold, 14)),
    TEXT(1.2*~coordinates(D), "D", FONT(Helvetica, bold, 14), COLOR(RGB, 0, 0, 1)),
    POLYGONS([coordinates(A), coordinates(B), coordinates(C)], COLOR(RGB, 0.9$3))
  ),
  plot3d([rhs~(eq__AD)[]], s=0..1, color=blue, thickness=5)
)

with_geom3d.mw

I'm not going to risk a stiff neck to read this but it sounds like a fixed point method.
Here is an example, adapt it to your own function f


 

Download fixed_point.mw

Look to the first line in the attached fie, this could also interest you

For lazzy person only

(PS: you will probably find everything you need in the Finance package; unfortunately, English not being my mother tongue and me not being at all familiar with financial terms, I used a summation to get the result. I have no doubt that if Finance has no secret for you you will find in this package the procedure that answers your question)

restart:
account_earns := 7.75/100;
years         := 30;
investment    := 4000;
                         0.07750000000
                               30
                              4000
add(Finance:-futurevalue(investment, account_earns, years-k), k=0..years-1)
                                      5
                        4.664152587 10 

But you must have been taught how to calculate this?
Does the term Geometric_series ring any bell?
 

Answer to your first question

This is a simple proposal (I used Maple 2015, newer versions produces more beautiful graphs).
I did not plot a ball at the vertices, but draw the path (sequence of edges) from the root of the tree to a termibal leaf).
A ball can be drawn (as you requested) but we need do know if you just want it to be displayed when it arrives on a vertex (simple) or if you want also bo follow it along the edges (more complex).
 

Download Animated_path_along_a_tree.mw


Maybe you could be interested in plotting different random paths?
Something like this....

f := proc()
  local U, Leaves, V, P, X, T2, dg, traj:
  U      := 1:
  Leaves := map(v -> if Degree(T, v)=1 then v end if, Vertices(T)[2..-1]);
  V      := combinat:-randperm(Leaves)[1];
  P      := ShortestPath(T, U, V());
  X      := GetVertexPositions(T)[P];
  traj   := CurveFitting:-BSplineCurve(X, t, order=1):
  traj   := subsop(3=(t=0..numelems(P)-1), traj):
  T2     := CopyGraph(T):
  HighlightVertex(T2, P, red); # highlight the vertices the ball hits
  dg     := DrawGraph(T2):
  display(dg, plot(traj, color=red, thickness=3, frames=numelems(P)));
end proc:

N := 10:
#animate(f, [], n=1..N, frames=N)

Once you know that your serie is the Taylor expansion of  1/(1+exp(x-5*t))^2 at t=0 (see here 232686-INFINITE-SERIES-SUM)
the answer is obvious

plot3d(1/(1+exp(x-5*t))^2, x=-10..10, t=-2..2, labels="x", "t", "u")

It's obvious that your serie can be written as a geometric serie (negative if lambda < 0, alternated if lambda > 0)

restart:
a := lambda*exp(t):
q := (-1)*lambda*(exp(t)-1):
U := (t, n) -> add(a*q^k, k=0..n):

#check
U(t, 4);

# let's set define S(t, lambda, n) this way
S := (t, lambda, n) -> lambda*exp(t) * add(((-1)*lambda*(exp(t)-1))^k, k=0..n):

# S(t, lambda, t, n) converges iif abs(q) < 1
# assuming lambda > 0 S converges iif abs(q) < 1 that is if 0 < lambda < 1/(exp(t)-1)
# Examples

T := 2.0:
plot([seq([k, S(T, 0.95/(exp(T)-1), k)], k=1..200)]);
plot([seq([k, S(T, 1.05/(exp(T)-1), k)], k=1..200)]);

So, what is lambda?

If the serie converges its limit is 

 a * 1/(1-q);
                         lambda exp(t)     
                    -----------------------
                    1 + lambda (exp(t) - 1)

This can be obtained by elementary calculus: the sum of a (converging)  geometric serie is  limit(a*(1-q^n)/(1-q) , n=+oo) = a/(1-q) for, because of the convergence condition abs(q)<1, q^n vanishes to 0.

after the line sys := ... write

indets(sys, name);
                        {a, b, c, d, x}

You have 4 equations and you claime having 3 unknowns.
Thus you miss one (a?) ... or maybe are you looking for another type of solution (the one that would minimize some functionnal of the 4 equations?

I took the liberty to consider a is an unknown (thus x is the only input).
A numeric solution can be obtained that way.

sys := unapply(sys, x):
# example for x=1
fsolve(sys(1))
    {a = -0.2722166549, b = -9.457425049, c = 2.635173636,  d = -3.098230243}

Is Maple able to find a formal solution of this same sys(1) system?

solve(sys(1), unknowns);

I waited about 5 minutes (Maple 2015) and couldn't get one.
So my guess is that Maple can't find a solution in terms of (here) x (alone)
 

Can you give the expression of the general term of this serie?
That is, if one write u(x, t)=u0(x) + u1(x)*t + ...un(x)*t^n+ ..., what is the expression of un(x)?
Or if you prefer, what is the recurrence relation between un(x) and un-1(x), un-2(x), ...un-p(x), (n>p)?

If we take the problem from the other side, one observe that the taylor expansion of u(x,t) with respect to t has its first terms equal to the first ones of the serie you give.

u := 1/(1+exp(x-5*t))^2:
S := map(simplify, taylor(u, t))

Of course it's easy do to what I'm going to do once we accept the previous result.

• the coefficients of the successive terms are 10, 50/2, 250/6, 1250/24, ...
  thus terms of the form 10 * 5^k / k! with k=0..+oo
   
• this suggest that the serie could be a Taylor expansion of some function w(x, t)
  moreover, giving the expression of the terms that contain t, this Taylor expansion, which should contain terms like (t-a)^k can be done at a=0
   
• for k=0 one has w(x, 0) = f(x) = 1/(1+exp(x))^2.
  let us observe a few derivatives of f(x) with respect to x:

  f := (1/(1+exp(x))^2:
  print~([seq(simplify(eval(diff(f, [x$k]), t=0)), k=0..5)]):
  

  we observe that these derivative are equal, up to a miltiplicative constant, to the terms that contain x in your serie
   
• if w(x, t) is of the form w(x+c*t), it's easy to see that 

  # set z=x+c*t
  diff(w(x+c*t), t$k) = diff(w(z), z$k)*c^k

  this could lead to the inference that your serie is the Taylor expansion (1/(1+exp(x+c*t))^2 around t=0
   
• under this hypothesis one easily finds that c=5

If the assumption that your serie is such a Taylor expansion, it converges for any real values of x and t.