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);  

Is this enough?

# p a polynom

type(p, polynom(anything, x)) and `not`(has(p, I))

returns true if p is a polynom in x and no coefficient is complex

I have often noticed that integration of a piecewise function can be difficult  (or even impossible), but that once the function is expressed using Heavisides,  the integration is done well.

uv := convert(u*v, Heaviside):
int(uv, x=-1..1, y=-1..1);
                               1 
                               --
                               12

I think it only be done by hand:

1. below the plot insert this

   P := [
   
   ];
   
   MyPoints := Matrix(nops(P)/2, 2, P);
   
   alignment := align=right:   # to be set in an ad hoc way
   plots:-display(
     p,
     seq(plots:-textplot(evalf[4]([MyPoints[n, 1], (MyPoints[n, 2])$2]), alignment), n=1..nops(P)/2),
     seq(plots:-pointplot([MyPoints[n, 1], MyPoints[n, 2]], symbol=cross, symbolsize=15, color=blue), n=1..nops(P)/2),
     seq(plot([[0, MyPoints[n, 2]], [MyPoints[n, 1], MyPoints[n, 2]]], linestyle=3, color=blue), n=1..nops(P)/2),
     seq(plot([[MyPoints[n, 1], 0], [MyPoints[n, 1], MyPoints[n, 2]]], linestyle=3, color=blue), n=1..nops(P)/2)
   ); 

2. right-click (control+click with Mac OSX) on the plot to select Prob Info / Nearest point on line
3. move the probe with the mouse to select a point
4. right-click (control+click) and select Probe info / Copy Data  (one single point at the times)
5. In the P:=[..] structure : paste the data (ctrl+v or command+v on Mac OSX)
6. go to step 3 to add a new point or to the next step
7. add a comma after each line but the last one in the P:=[..] structure
8. replace the decimal separator (comma) by a point
9. execute the block 

probe.mw
 

What is your definition of "minimal set"?

If you mean "the set which contains the minimum number of sets", then your second example produces a minimal set among others.

For instance the set 

{{}, {5}, {1, 6}, {4, 6}, {2, 3, 6}}

has also 5 elements whose unions generate all intersections of the elements the set of your second example

{{1, 2, 3, 6}, {1, 4, 5, 6}, {2, 3, 4, 6}}

contains.
 

Here is an idea to obtain "a set which contains the minimum number of sets".

S := {{1, 2, 3}, {1, 4, 5}, {2, 3, 4}};
N := nops(S);

eS   := op~(S);
T__0 := {{}, (`{}`~(eS))[]};

cross := table([]):
for m from 1 to N-1 do
  for n from m+1 to N do
    mn := [m, n]:
    cross[mn] := S[m] intersect S[n];
    print(mn, %);
  end do:
end do:

T__1 := {{}, entries(cross, nolist)};
T__2 := T__1 union {(eS minus op~(T__1))};

minimal_sets.mw

REMARK 1
Replacing the line

T__1 union {(eS minus op~(T__1))}

by

T__1 union {(eS minus op~(T__1))} minus {{}}

gives an "even more minimal set" which has the good properties... at least if no intersection between two elements of the initial set is the empty set.


REMARK 2
My method doesn't return a correct answer if the initial set contains an element U whose intersection with any other element is empty.
But improvements could be done

As @dharr mentioned it, using useassumptions is extremely slow (I stopped the computation even before it completed).

Here is a workaround.

STEP 1: discard obviously wrong solutions

nops(sol4cons):
                               18

# A simple way to discard solutions with b=c
sol4cons_cleaned := select(`not`@has, sol4cons, b=c):
nops(%);
                               17

# Unfortunately this works only if the solution contains the relation b=c.
# It doesn't discard sol4cons[6] which contains b=0 and c=0.
#
# A more general method
sol4cons_cleaned := map(s -> if evalb(eval(b<>c, s)) then s end if, sol4cons):
nops(%);

# numbers discarded are 

sol4cons_cleaned := map(n -> if `not`(evalb(eval(b<>c, sol4cons[n]))) then n end if, [$1..nops(sol4cons)]);
                               16
                             [4, 6]

The problem is that the remaining solutions are not correct, for instance 

sol4cons[1]
 /                              1           2                  
{ a = a, b = b, c = c, lambda = - alpha A[0] , mu = mu, p = p, 
 \                              3                              

                                  \ 
  A[-1] = 0, A[0] = A[0], A[1] = 0 }
                                  / 

and nothing insures that a^2-b^2+c^2 > 0.


STEP 2: do a second pass on sol4cons_cleaned
So there is another workaround:
For each remaining solution in sol4cons_cleaned : "re-solve" the solutions with the constraints a^2-b^2+c^2 > 0 and b<>c:

goodsols := map(s -> solve(s union {b<>c, a^2-b^2-c^2>0}), sol4cons_cleaned):

There remain two problems

• Two  solve fail:

• n := 14: # and n:=15
  sol4cons_cleaned[n]:
  solve(sol4cons_cleaned[n] union {b<>c, a^2-b^2-c^2>0})
  Error, (in unknown) numeric exception: division by zero
  

I guess this mean that "preliminary" solutions number 14 and 15 cannot be solutions that respect the constraints, but I'm not at all sure of this interpretation
 

• the 5th solution takes a long time to be re-solved (I stopped the process)
   

Final remark
I guess useassumption leads to extremely time consuming computations because there is a huge number of solutions to compute.
For instance, re-solving sol4cons_cleaned[1] generates 12 solutions:

n := 1:
sol4cons_cleaned[n]:
[solve(sol4cons_cleaned[n] union {b<>c, a^2-b^2-c^2>0})]:
nops(%);
                           12
print~(%%)