@Michael_Watson 

Download asympt.mw

One way with geom3d package

EDITED VERSION

 

restart:
with(geom3d):

# first line

plane(P1, 16*x-2*y-11*z, [x,y,z]):
plane(P2, 14*x-y-10*z-3, [x,y,z]):
line(L1, [P1,P2]):
Param_Eq__1 := Equation(L1,'t');
                    [1                     ]
                    [- + 9 t, 4 + 6 t, 12 t]
                    [2                     ]
# second line (=L2)
# let u be the common value of (x-2)/3, (y-5)/2 and (z-2)/4

X := solve((x-2)/3=u, x);
Y := solve((y-5)/2=u, y);
Z := solve((z-2)/4=u, z);
r := u=solve(X=rhs(Param_Eq__1[1]), u);
eval([X, Y, Z], r)
                               1      
                         u = - - + 3 t
                               2      
              [1                     ]
              [- + 9 t, 4 + 6 t, 12 t]
              [2                     ]

Thus lines L1 and L2 are  identical.

First problem only (I'm going to bed)

Step 1: define M and N this way

restart:
M := (a, b) -> (a+b)/2;
N := (a, b) -> a*b/M(a, b);

Step 2: the recurrence is

a(n+1) = M(a(n), b(n)),  b(n+1) = N(a(n), b(n))

Set  z(n+1) = a(n+1) * b(n+1) and observe that z is invariant under the transformations M and N.

z(n+1) = a(n+1) *  b(n+1) = a(n) * b(n) = ... = a(0) * b(0)

Step 3 : find a limit (aka stationary) solution by solving this system of equations (could be done by hand)

solve({M(x,y)=x, N(x, y)=y})
                         {x = y, y = y}

Thus a stationary solution is such that x=y (or limit(a(n), n=+oo) =  limit(b(n), n=+oo) if you prefer


Step 4 : inject this stationary solution within z : 
As limit(a(n), n=+oo) =  limit(b(n), n=+oo), one gets limit(z(n), n=+oo) = z(0) = a(0)*b(0) and thus 

limit(a(n), n=+infinity) = limit(b(n), n=+infinity) = sqrt(a(0)*b(0))

Example : a(0)=1 and b(0)=2 gives  limit(a(n), n=+oo) =  limit(b(n), n=+oo) = sqrt(2)

Does this suit you?

(u instead of u(y))

restart
alias(u=u(y, z));
de := diff(u, y$2)+3*delta*diff(u, y)^2*diff(u, y$2)=p;

# and for subsequent computations
eval(de, u=U(y));
# example
dsolve(%, U(y));

Note that the trick alias(u(y)=u(y, z)) won't produce the desired result

It seems odd to want to create a diagonal system but you probably have your own reasons.
Here is a possibility

restart
interface(version)
Standard Worksheet Interface, Maple 2015.2 ....

with(LinearAlgebra):

A := (z, M) -> DiagonalMatrix([z$M]):
V := (n, M) -> Vector(M, i -> u[i, n]):

# example of use,
# the neutral operator &. to avoid the evaluation of the matrix vector product


A(omega[2]*(2-b[1]), 3) &. V(2, 3);

# your question:
N := 4:
m := 3:

for j from 2 while j <= N do
   printf("%s j = %d %s\n", cat("-"$50), j, cat("-"$50));
   A(omega[2]*(2-b[1]), m) &. V(j, m)
   +
   A(2*b[1]*omega[2]-b[2]*omega[2]-omega[2]+1, m) &. V(j-1, m)
   -
   (sum(A(omega[2]*(b[l+2]-2*b[l+1]+b[l]), m) &. V(j-l-1, m), l = 1 .. j-2))
end do;

REVISED 2021/8/06
A revised version inspired by  the answer acer made to this question 232636-Is-It-Possible-To-Construct-Binary-Operators

System.mw
 

Your parameterization seems to be too complex for Maple to be able to handle the assumptions.
I suggest you to use this one

assume(s < t, t < 0, 0 < u, u < v);
r := x -> piecewise(
            s <= x and x <= t, b*(1 -((t-x)/(t-s))^(2)),  
            t <  x and x <  u, b,  
            u <= x and x <= v, b*(1- ((x-u)/(v-u))^(2))
          )

Where s, t, u and v correspond to

• s -->  -a
• t -->  -a + epsilon
• u --> +a - varepsilon
• v --> +a

r(s);
                               0
r((s+t)/2):        
r(t);
                               b
r(0);
                               b
r(u);
                               b
r((u+v)/2):
r(v);
                               0

See in the attached file
reparameterization.mw


PS: I thought you caps where ellipses with half small axes equal to epsilon or varepsilon (depending if they cap the left or the right section of the cylinder) and half great axes equal to b?
This not the case xith what you did (the tangents to the elliptic cap at x=s or x=v are not vertical.
Maybe you wanted this?

r := x -> piecewise(
            s <= x and x <= t, b*sqrt(1 -((t-x)/(t-s))^(2)),  
            t <  x and x <  u, b,  
            u <= x and x <= v, b*sqrt(1- ((x-u)/(v-u))^(2))
          )

LAST POINT: since you mentioned cylinder, I don't know if you are aware of the SurfaceOfRevolution  function?
If not you might like this

Student:-Calculus1:-SurfaceOfRevolution(eval(r(x), [s=-3, t=-2, u=2, v=3.5, b=2]), x=-3..3.5, output=plot, scaling=constrained);

To avoid the errors just change v__2p := 0.0 into v2p := 0*Unit(('m')/('s'));

zadanie_z_jednostakim_-_problem__rep.mw

I think you might also like v__1k and v_2k  to be expressed with their units (m/s)?
In such a case here is a solution

m__2 := .4*Unit('kg');
m__1 := .3*Unit('kg');
x__w := .7*Unit('m');
v__2p := 0.*Unit('m'/'s');
v__1p := 2*Unit('m'/'s');

pot := m__1*v__1p+m__2*v__2p = m__1*v__1k+m__2*v__2k;
kin := (1/2)*m__1*v__1p^2+(1/2)*m__2*v__2p^2 = (1/2)*m__1*v__1k^2+(1/2)*m__1*v__2k^2;

solve([pot, kin], [v__1k, v__2k])

Here is a solution: at the end of uor code you can include this (see the attached file for a complete worksheet)

start := StringTools:-Ord("A");
letters := StringTools:-Char~([$start..start+25])

_iquo := iquo(n, 26):
_irem := irem(n, 26):

if n <= 26 then 
  LETTERS := letters[1..n]
else
  LETTERS := letters[]:
  for i from 2 to _iquo do
    LETTERS := LETTERS, cat~(letters, i-1)[]
  end do:
  LETTERS := LETTERS, cat~(letters[1.._irem], _iquo)[]
end if:

newlist       := [LETTERS] =~ z;            # elements of the form "A2"=[1, 2]
other_newlist := parse~([LETTERS]) =~ z     # elements of the form A2 =[1, 2]

letters.mw


REMARK : there is no real assignment here in the sense of ?assign.
The reason is that some "letters" such as D, I, O are protected, so you can't assign them unless unless you define them as local  (which could be dangerous).
An example is provided at the end of the attached file

Simple reasoning:
I assume that all the points have integer coordinates (right?) and that mov1ng from A to B means "always move in the same horizontal direction and in the same vertical direction" (if it's not the case a pathh from A to B can have an infinite length).
More of this I understand that you want to go from A to E by going fromA to B, next from B to E and so on (right?).
If I'm wrong on this last point a minor modification of what I'm going to say can be derived to find a shorter path that visit all the points once/

Let (x_A, y_A) the coordinates of A and (x_B, y_B) those of B.
Let  x_m=min(x_A, x_B), x_M=max(x_A, x_B), y_m=min(y_A, y_B), y_M=max(y_A, y_B).
Let h_x = x_M-x_m and h_y = y_M-y_m.
Note that h_x  and h_y  are 2 integers; thus there exist exactly  C(h_x+h_y, h_y) = (h_x+h_y)!/(h_x!_*h_y!) paths from A to B that have all the same length h_x + h_y .
Once arrived at B, do exactly the same thing to determine the number of paths from B to C and their length
(which is just abs(x_C-x_B)+abs(y_C-y_B))

Thus the length of any of all the paths from A to E is equal to 
abs(x_B-x_A)+abs(y_B-y_A) + .... + abs(x_E-x_D)+abs(y_E-y_D) 
There exist exactly C(abs(x_B-x_A)+abs(y_B-y_A), abs(y_B-y_A)) _* ... _* C(abs(x_E-x_D)+abs(y_E-y_D), abs(y_E-y_D)) such paths

For instance, if A=(0,0) and B = (2, 3), there exist 6 paths ((2+3)!/2!/3!) of length 5 (X denote a move of length in the horizontal direction and Y a move of length 1 in the vertical direction)
X X Y Y Y
X Y X Y Y 
X Y Y X Y
X Y Y Y X
Y X X Y Y 
Y X Y X Y
Y X Y Y X 
Y Y X X Y
Y Y X Y X
Y Y Y X X


For the graph you give there are:

• 2 paths of length 2 from A to B
• 5 paths of length 5 from B to C
• 10 paths of lrngth 5 from C to D
• 2 paths of length 2 from D to E

and thus 200 shortest paths of length 14 of the form A->B->C->D->E

Note that the paths of type A->B->D->E->C have length 13 (see above my second assumption).

In case you want to find a cycle which passes from all the points (not a path), which means the ending point is also the starting point, it's bovious to see that the shortest cycle is A->B->D->E->C->A (3360 such cycles if we do not account for the starting point) with length 20

PS:the distance you refer to is the Mahalanobis distance, also dubbed the "taxi distance" by reference to the path a taxi takes in New York downton where streets are either pairwise parallel or orthogonal.

use combine(expression, 'units')

2*Unit('mA') * 3*Unit('V');
              6 Units:-Unit('mA') Units:-Unit('V')

combine(%, 'units');
                       3                  
                      --- Units:-Unit('W')
                      500                 

and thus 

local I:

I := Vector(3, i ->(2*i)*Unit('mA'));
U := Vector(3, i ->(3*i)*Unit('V'));
P := combine(I*~U, 'units');

Your ode seems strange. In terms of if..then.. else structure it means (if I'm not mistaken)

if abs(x(t)) > 0 then 
   something
elif abs(x(t)) < 1 then 
   something else
end if

What happens if abs(x(t)) > 1?
Are you sure the second condition is not 

elif abs(x(t)) > 1 then 
   something else
end if

???

In order to show you how we can handle this type of problem I decided arbitrarily to solve a problem where:

• when x(t) = -1 or +1 then the source term changes to 0
• when x(t) = 0 then the source term changes to 2

What I would do (there are probably other possibilities) is to introduce a discrete_variable A(t) (the source term is 2*A(t) where A(t) is either 0 or 1 depending on the value of x(t))

restart:
ode := diff(x(t), t)=v(t), diff(v(t), t) - 2*A(t) = 0;

# set some initial conditions, for instance

ic  := x(0)=-2, v(0)=1, A(0)=1;               
sys := {ode, ic};


sol := dsolve(
  sys, 
  numeric, 
  discrete_variables=[A(t)], 
  events=[[x(t)=1, A(t)=0], [x(t)=-1, A(t)=0],[x(t)=0, A(t)=1]]
):
proc(x_rkf45)  ...  end;

plots:-odeplot(sol, [t, x(t)], t=0..2);
plots:-odeplot(sol, [x(t), v(t)], t=0..2)

Download discrete_variable.mw

If I didn't correctly understand your problem, forget this reply.

The point you search is the "double point" of the curve.
That is in fact 2 points M1=(f(t), g(t) and M2=(f(tau), g(tau) which verify these three relations

{f(t)=f(tau) g(t)=g(tau), t <>tau}

For unicursal curves the method to find these points relies upon a reparameterization  in terms of P and S where

P=t+tau
S=t*tau

Your curve is of this type (Unicursal_curve).
The first step consists in transforming the system 

{f(t)=f(tau) g(t)=g(tau), t <>tau}

into a system 

{F(S, P), G(S, P)}

In the second step this system is solve (analytocally on some cases, but numerically in yours) in P and S.
In the last step one finds t while solving the equation

t^2 - S*t + P = 0

In your case this gives

f := t -> exp(t) + exp(-t^2 - 2*t):
g := t -> exp(-t^2) + exp(-t^2 - 2*t):

e1 := simplify(expand(f(t)-f(tau))/(t-tau));
e2 := simplify(expand(g(t)-g(tau))/(t-tau));

parameterization := [t, tau] =~ [solve(t^2-S*t+P, t)]:

newsys := eval([e1, e2], parameterization):
PS_sol := fsolve(newsys);
             {P = -0.1824567273, S = -1.566488415}

t_sol := fsolve(eval(t^2-S*t+P, PS_sol), t)
                   -1.675392298, 0.1089038833

double_point := [x,y] =~ [f,g](t_sol[1]);  # same values for tsol[2]
               [x = 1.909852743, y = 1.783007571]

double_point.mw


Remark: this method considers turning points as  double points, which means (it's not your case) that one mist verify that the points found are real double points.

Forget it, the expression is far too complicated for Maple (or any other CAS) to ever find its inverse Fourier transform.
Personnaly I would try a discrete approach and do something like this

# evaluate u_tilde at x=L

uL_tilde := eval(u_tilde, x=L):

# verify that the only indeterminate in uL_tilde is omega

indets(uL_tilde, name);


# transform uL_tilde as a function of omega

uL_tilde := unapply(uL_tilde, omega):

                            {omega}
with(DiscreteTransforms):
A := Vector(1024, [seq(evalf(eval(uL_tilde, omega=n*30/1024)), n=1..1024)], datatype=complex[8]):
B := InverseFourierTransform(A):

plots:-display(
  plot([ seq([t-1, Re(B[t])], t=1..1024) ], color=blue, legend="real"),
  plot([ seq([t-1, Im(B[t])], t=1..1024) ], color=red , legend="Img")
);

of course there remains many thinfs to fix:

• The omega range (implicitely set to 0..30 in the example) is to be set to something more physical.
• The number of discretization points (1024) should be adjusted too.
• The final plot is done with a wrong absissa (I took the values 0, 1, ... 1023 as they hsould be of the form 0, C, 2C, ... 1023C, where C depends on the range of omega and the definition of InverseFourierTransform [there must be some Pi somewhere and probably a square root too])
   

I don't know if this will help you.

Also take a look to the package SignalProcessing, it could be useful.

A few remarks:

• If what matters are the coefficients of the fitted model the simplest way to get them is

  Para_Weight_Coeff := NonlinearFit(Para_func, map(convert, Weight_time, unit_free), map(convert, Pull_Weight, unit_free), x, output = parametervalues)
  
              [A_Const = 1.8619791553713554, C_Const = 6.034041446128156]
  

• For the reference vocabulary see "International vocabulary of metrology – Basic and general concepts and associated terms (VIM)", f0e1ad45-d337-bbeb-53a6-15fe649d0ff1

  
  The error you got and which forces you to some gymnastics is that your model is not physical from the point of view of dimensional analysis.
  To explaint that, take the logarithm of both sides of your model to make it linear with respect to x  (provided A_Const, C_Const and x are both positive):

  ln(Para_func) = ln(A_Const)+x*ln(C_Const)

  From dimensionnal analysis, A_Const and  Para_func (in fact Pull_Weight)  must have the same dimension.
  The case of the term  x*ln(C_Const) is more complex for it must be dimensionally homogeneous to ln(Para_func). Phisically this have sense only if x=1 and C_Const and  Para_func both have the same dimension, or if x and  ln(Para_func) both have the same dimension and C_Const is a number.
  To understand where the error comes from just to do this

  LinearFit(ln(A_Const)+x*ln(C_Const), Weight_time, `~`[log](Pull_Weight), x)

  You will obtain the same error, and again the same error if you just do this

  ln(size1);
  Error, (in Units:-Standard:-ln) the function `Units:-Standard:-ln` cannot handle the unit m with dimension length
  
  
  cos(size1);
  Error, (in Units:-Standard:-cos) the function `Units:-Standard:-cos` cannot handle the unit m with dimension length
  
  
  (size1)^(sqrt(2));
  Error, (in Units:-Standard:-^) a unit can only be raised to a rational power
  

  This is simply because dimensions always combine multiplicatively through terms like 

  Q__1^alpha__1 * Q__2^alpha__2 *...
  

  where Q__n are quantities (length mass) whose expressions use Units, and thankfully Maple knows that
  Thus the error you got has nothing to do with NonlinearFit: it's just physical inconsistency that Maple warns you about.
  
  Of course you may have quantities within a mathematical functions like exp, ln, cos and so on, but these quantities can only appear through multiplicative non dimensional expressions like Q__1^alpha__1 * Q__2^alpha__2 *... 
   
• When you write 

  Eq1 := Para_Weight_Coeff1*Unit('kg'/'s')*Para_Weight_Coeff2^t*Unit('s')

  this means that Para_Weight_Coeff2^t is a quantity homogeneous to a time, but later you say that t is itself a time, which is a nonsense.
  You may give the feeling that Eq1 is a quantity homogeneous to a mass (happily Pull_Weight is a mass too ;-) ) but there is something fundamentally wrong in this construction.

Here is a detailed answer which took me a lot of time to write (more or less) correctly.
I met two problems:

1. first, your demand is not clear (not to say the transition matrix you use is unnecessary complicated),
2. secondly I don't know what your knowledge about stochastic processes is.

Please look to the content of the attached file.
It contains many things that people familiar with stochastic processes should understand. Some of these things are objective, while others are purely subjective and reflect my own understanding of your often too hazy questions.

Markov_process.mw