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

EDITED (see WARNING below)


As you have probably noticed the reason is that fsolve returns a set while Isolate returns a list.
A set with elements such that a=b is ordered according to the left hand side terms... and after x1 comes x10, and x11, ...x19, and then only x2.
This disturbs me too and I use to use this trick to reorder the output of fsolve in a more "natural" (?) order

fsol := fsolve([f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12]); 

fsol_order     := map(u -> parse(substring(u, 2..-1)), lhs~([fsol[]]));
natural_order  := sort(fsol_order, output=permutation);
natural_vars   := lhs~([fsol[]])[natural_order];
natural_fsolve := natural_vars =~ eval(natural_vars, fsol);

WARNING
as suggested by @acer , I have to say that there are more concise and efficient ways to do this.
Maybe I misunderstood what exactly you were looking for? All the more so that @tomleslie's answer seems to satisfy you and that acer considers it the right answer.
What I had understood is that you wanted fsolve to return a result "in the same order" as Isolate? Or at least to rearrange the output of fsolve in order it lokks like the output of Isolate?
A simpler method than the one above to do this is

# define some "order" (here the order you used in Isolate)
vars := [x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12]:

ans2 := RootFinding[Isolate]([f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12], vars)[4]

ans1 := fsolve({f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12}); 
# reorder the outputs of solve given the "order" vars
vars =~ eval(vars, ans1)

But I may be off the mark.

A lot of syntax errors (corrected in the attached file).
The solutions contain a singularity around T=0.0063666816, so I ploted them on a restricted range compared to yours (T=0..12)

Watchout : some parameters are not present in the differential system: dsolve can handle this but it would be better to suppress them.
Use this command

indets({ODEs, bcs}, name) minus{T}

to see what are the parameters of the system.

Here is an excerpt of the attached file

dSol:= dsolve({ODEs, bcs}, numeric, parameters = lhs~(params)):

Warning, the unknowns {__gamma, alpha, beta__o, psi, sigma, xi} were specified in 'parameters', but are not present in the input system


getPlot := proc( pn::name, pv::list, fv::function, Trange )
  local where_is_it, MyColor, plts, val, test_case:
  where_is_it := ListTools:-Search(pn, lhs~(params)):
  MyColor := () -> ColorTools:-Color([rand()/10^12, rand()/10^12, rand()/10^12]):
  plts := NULL:
  for val in pv do
    test_case  := subsop(where_is_it = val, rhs~(params)):
    dSol(parameters=test_case):
    plts := plts, plots:-odeplot(dSol, [T, eval(fv)], T = Trange, color=MyColor(), legend=typeset(pn=val) ):
  end do:
  return plots:-display(plts):
end proc:

getPlot(M__h, [0.152, 0.252, 0.352, 0.452, 0.552], B(T), 0..0.006 );

Errorplot_corrected.mw


Here is an updated version of the Errorplot_corrected.mw which enables to plot the solutions for combinations of several parameters.

Errorplot_updated.mw

The idea is to print the procedure 'PlanePlot' into an mpl file using writeto:
(the code below comes from acer's answer to one of my questions, see here 221000-How-To-Redirect-Showstat-Into-A-File-)

restart:
with(Student[LinearAlgebra]):
interface('prettyprint'=1):
interface('verboseproc'=3):
f := cat(currentdir(), "/Desktop/test/PlanePlot.mpl");
writeto(f);
print('PlanePlot');
writeto('terminal'):

Next, click on Open in the main menu bar and open the mpl file  PlanePlot.mpl.
The procedures appears in a new worksheet.
Insert 

restart:
with(Student[LinearAlgebra]):
MyPlanePlot := 

before the name proc.
In Maple 2015 I have to do some slight modifications to make the code syntactically correct.
Basically I replaced thins like

B := map(proc(t)
option operator, arrow;
      :-LinearAlgebra:-Multiply(t, (max(1.0, 
      :-LinearAlgebra:-Norm(p, 2, 'conjugate' = false)))/(
      :-LinearAlgebra:-Norm(t, 2, 'conjugate' = false)))
    end proc;, B);

by things like

B := map(t -> 
      :-LinearAlgebra:-Multiply(t, (max(1.0, 
      :-LinearAlgebra:-Norm(p, 2, 'conjugate' = false)))/(
      :-LinearAlgebra:-Norm(t, 2, 'conjugate' = false)))
      , B);

Insert DEBUG(): after the declaration of local variables.
Here is the final result
PlanePlot.mw

Run the command

MyPlanePlot( -3*x+2*y+z = -3, [x,y,z], normaloptions=[shape=harpoon] )

(for instance).
You then enter the debug mode and can follow step by step what PlanePlot does.

For this example the main operations are
 

P0:=lhs(P) - rhs(P);

P0,d:=selectremove(has,P0,vars); 
d:=evalf(-d);

p:=d*n/:-LinearAlgebra:-Norm(n,2,'conjugate'=false)^2; 
p0:=p; 
B:=:-LinearAlgebra:-NullSpace(Matrix([n[1],n[2],n[3]],'datatype'=float[8])):
B:=map(t->:-LinearAlgebra:-Multiply(t,max(1.0,:-LinearAlgebra:-Norm(p,2,'conjugate'=false))/:-LinearAlgebra:-Norm(t,2,'conjugate'=false)),B)

Changing numpoints is a way, another one is to use the option adaptive=true.

Look at the attached file to understand why the graphics are not the same and see how to solve the problem.

plots.mw

Here is the way which mimics what you say (even if the things could be done differently)
The core is this procedure

operations := proc(A)
local V := `<|>`([A + 2, 2*A, A^2, A mod 2, A + 3, 3*A, A^3 ,A mod 3]):
if _npassed = 1 then
  return V
else
  return < _passed[2], V >
end if:
end proc:

It could (should?) be modified to use a type checking of the second argument when it's passed.
See how to use this procedure in the attached file and how to export the final result.

operations.mw

One uses BarGraph :

BarChart(dataset, legend = ["vec 1", "vec 2"], axis[2]=[tickmarks=[seq(i-1/2=Variables[i], i=1..numelems([Variables]))]]);

but the labels are written on the bottom of the boxed axes and I don('t know how to place them on the y-axis???

the second uses an ad hoc way which seems more flexible

plots:-display(
  seq(plottools:-rectangle([0, n-1/3], [dataset[1][n], n], color=green), n=1..N),
  seq(plottools:-rectangle([0, n], [dataset[2][n], n+1/3], color=gray), n=1..N),
  axis[2]=[tickmarks=[seq(i=Variables[i], i=1..N)]]
)

You will find a few customizations in the attached file plus a transposition to  ColumnGraph
barchart.mw

A solution.

EDITED 6/9 9:39 am

Transform the fsolve  problem (which is the one who fails) into a minimization one
 

Maple 2015.2 (iMac OSX Mojave)

b       := 0.1e-3:
epsilon := 0.1e-4:
theta   := .5:
n       := 1000:
alpha   := .1:

CodeTools:-Usage( f(b, epsilon, theta, n, alpha) )
memory used=1.27GiB, alloc change=134.01MiB, cpu time=4.38s, real time=3.39s, gc time=1.75s
                0.831910275472156306816487652377
alpha   := .01:

CodeTools:-Usage( f(b, epsilon, theta, n, alpha) )
memory used=69.18MiB, alloc change=8.00MiB, cpu time=2.10s, real time=2.04s, gc time=109.54ms
                0.981860839533060088397020565756
alpha   := .001:

CodeTools:-Usage( f(b, epsilon, theta, n, alpha) )
memory used=83.73MiB, alloc change=256.00MiB, cpu time=2.22s, real time=2.15s, gc time=111.72ms
                0.998177604561510290591588684447
alpha   := .0001:

CodeTools:-Usage( f(b, epsilon, theta, n, alpha) )
memory used=86.21MiB, alloc change=0 bytes, cpu time=2.26s, real time=2.18s, gc time=185.59ms
                0.999817676230524565857083211528

Here is the code
Minimize.mw

A few remarks :

1. In the first par of the attached file I kept using int(PDF(RV, x), x=0..a) but it's more clever to replace these by eval(CDF(RV, x), x=a) : indeed, the expression of the CDF is (generally) a part of the definition of the random variable.
   This means integration is not needed ; use this to verify that the cdf is not obtained by integration of the pdf (or the latter by derivation of the former):
    

   B   := RandomVariable('Beta'(a, b)):
   L   := [attributes(B)][3];
   exports(L);
   L:-CDF(x)
   

   
   With this trick the computational is about 10 times faster:
   (note that the results can even been obtained with Digits:=10)
    

   b       := 0.1e-3:
   epsilon := 0.1e-4:
   theta   := .5:
   n       := 1000:
   alpha   := .1:
   
   CodeTools:-Usage( g(b, epsilon, theta, n, alpha) )
   memory used=17.02MiB, alloc change=0 bytes, cpu time=134.00ms, real time=134.00ms, gc time=0ns
                   0.831910275472156306816487652377
   alpha   := .01:
   
   CodeTools:-Usage( g(b, epsilon, theta, n, alpha) )
   memory used=22.36MiB, alloc change=0 bytes, cpu time=170.00ms, real time=214.00ms, gc time=0ns
                   0.981860839533060088397020565756
   alpha   := .001:
   
   CodeTools:-Usage( g(b, epsilon, theta, n, alpha) )
   memory used=35.32MiB, alloc change=0 bytes, cpu time=421.00ms, real time=327.00ms, gc time=149.80ms
                   0.998177604561510290591588684448
   alpha   := .0001:
   
   CodeTools:-Usage( g(b, epsilon, theta, n, alpha) )
   memory used=37.16MiB, alloc change=0 bytes, cpu time=273.00ms, real time=291.00ms, gc time=0ns
                   0.999817676230524565857083211528
   

2. I wrote 'Beta' instead of the lengthy BetaDistribution
    
3. More important: the assuming hypothesis

   phi := PDF('Beta'(alpha, beta), x)  assuming x >= 0, x < 1:

   in the body of procedure f are important because the PDF is always a piecewise function for a random variable with bounded support... and integrations of piecewise functions is generally costly.
   For the same reason I used 

   r := CDF('Beta'(alpha, beta), x)  assuming x >= 0, x < 1: 

   in the body of procedure g, to ease Optimization:-Minimize.
    
4. The two last commands in procedure g are equivalent to

   B := RandomVariable('Beta'(alpha, s+n)):
   return Probability(B <= b, numeric)
   

    

A way

Step one: use writeto to redirect the output in some text file
Step two: read the file and extract the desired information
 

restart:
with(Student[LinearAlgebra]):
infolevel[Student[LinearAlgebra]] := 1:

MyFile := cat(currentdir(), "/PP.txt"):
writeto(MyFile):

PlanePlot( -3*x+2*y+z = -3, [x,y,z], normaloptions=[shape=harpoon] ):
close(MyFile)


writeto(terminal)
with(StringTools):

line := readline(MyFile):
while line <> 0 do
   last := line;
   line := readline(MyFile);
   if line <> 0 and Search("basis vectors:", line) =1 then
     bv := StringSplit(line, "basis vectors:")[2];
     bv1, bv2 := parse(%);
   end if:
end do:

`basis vector 1` = bv1; 
`basis vector 2` = bv2;

PlanePlot.mw

tpn(...) isn't defined.

See here that it should work once tpn(...) is defined.
tpn.mw

Only the first one solved, just do the same thing for the others/

Remarks

1. dL/dt in Maple is not at all the derivative of L(t) with respect to t.
   Look at ?diff
2. Are T and R constants or functions of t?

You will find below a correct syntax in both cases : R = Cste and R = R(t)

If you are lazzy to write R(t) or if you just want a simpler expression you can use  aliases (?alias)         
 

restart:
alias(L=L(t)):
eq__1 := diff(L, t) = -L*b[2]+R*g[1];

sol__1 := dsolve(eq__1);

#-----------------------------------------------------------
      
restart
alias(L=L(t)):
alias(R=R(t)):

eq__2 := diff(L, t) = -L*b[2]+R*g[1];

sol__2 := dsolve(eq__2)
                
# 2 examples

eval(sol__2, [R=t, Int=int]);
eval(sol__2, [R=cos(t), Int=int]);

Other tips (?)

alias(`dL/dt`= diff(L(t), t)):
ode := `dL/dt`= -L*b[2]+R*g[1];
                    dL/dt = -L b[2] + R g[1]

... and if you REALLY want to write the equations as you did, you can do this:

restart:
alias(L=L(t)):
alias( `#mfrac(mo("dL"),mo("dt"))` = diff(L(t), t) ):

eq := `#mfrac(mo("dL"),mo("dt"))` = -L*b[2]+R*g[1]; #look at the output

# next
dsolve(eq);