A plottools object (like disk as proposed by @dharr) is always behind a plot object whatever the order you use in display.

A workaround is to use only plot objects, for instance :

p0 := plot([[x0, y0]], style=point, symbol=solidcircle, symbolsize=100, color=black);
p1 := plot([[x3, y3]], style=point, symbol=solidcircle, symbolsize=100, color=black);
display(c,p0,p1,size=[300,300]);  # place "c" behind the two big black points

Another way which still uses ploottols:-disk consists in drawing the blue arc only toutside the disks, let's say from t=first to t=last, where tfirst > 0 and tlast < 1.
Here is the corresponding code.
 

first := select(is, [solve((f(t)^2+g(t)^2)=0.02^2)], positive)[]:
last  := 1-first:
c     := plot([f(t),g(t),t=first..last],thickness=2,scaling=constrained, color=blue,axes=none):
p0    := disk([x0,y0],0.02,color=black):
p1    := disk([x3,y3],0.02,color=black):
display(c, p0, p1);  # same rendering than display(c, p0, p1);

Note that the plotting order within display is now irrrelevant.

Disk+Bezier.mw

There is no such command, and it's a good thing.

My feeling while saying this is that experimental data may have very different forms (a continuous analogic record, a digital sample), can be error free or not, the comparison criteria are generally proble-dependent and user-dependent, ...
More of this what are you going to do with a value of the discrepancy between experiment and simulation ? Align or Calibrate the simulation, or just do nothing more ?
So, IMO, it's better than Maple doesn't propose anything because it's likely that no one would be satisfied.

To summarize, you have todo the job yourself.
To help you doing this, here is a notional example where:

• the experimental data are measurement errors free,
• the experiment perfectly represents the reality (as a rule of thumb this is never the case)
• the mathematical model used in the simulation differs from the reality (that is from the model which mimics the reality perfectly).

One main poin I emphasize in this worksheet, is that the choice of a relevant "metrics" to characterize the discrepancy between experiment and simulation is up to the user.

In this worksheet I use two free parameters in the model (g and H) that could be adjusted to lessen the discrepancy between the experimental data and the simulation results.

The las criterion is base on the chronometry discrepancy, for his if often something which is of great concern (from my experience). Once you know how to compute it (procedure find(h)), you can imagine any criterion you want based on this discrepancy.

Download Experiment_Simulation_Comparison.mw

I use to use Maplets and we must be very carreful when doing this.
From my experience:

• I recommend to ALWAYS load the package Maplets
• I advice you to use with at the highest level and avoid uses within procedures

Here is your file corrected (the print where used to see what Display was responsible for the error [this was the last one])
MapletsDisplayError_Corrected.mw

PS: I use Maple 2015.2 and I hope this file still works for newer versions.
       (of course I answered NO two times to the checkForFile maplets for I don't have any of the
      requested files)

As Tom sais, the message is pretty clear. When this happens the good  thing to do is this

indets(d1, name); # for instance

To fix this, replace the last line 

plot([d1, d2, d3], eta = 0 .. 1)

by

plot(eval([d1, d2, d3], j), eta = 0 .. 1)

I guess that minimize has some difficulties in handling infinite search domains (which should exclude using [x > 0, y> 0, z > 0]) 
Look at this

restart;
sol := M -> minimize(x^3 + y^3 + z^3 - 4*x*y*z, [x=0..M, y=0..M, z=0..M], location);
                    / 3              3    3  
       M -> minimize\x  - 4 x y z + y  + z , 

                                                       \
         [x = 0 .. M, y = 0 .. M, z = 0 .. M], location/
sol(10.^100)
                      300   /[ /                       100  
  -1.00000000000000 10   , { [{ x = 1.00000000000000 10   , 
                            \[ \                            

                           100                         100\   
    y = 1.00000000000000 10   , z = 1.00000000000000 10    }, 
                                                          /   

                        300]\ 
    -1.00000000000000 10   ] }
                           ]/ 
sol(+infinity)
                0, {[{x = 0, y = 0, z = 0}, 0]}

In your code the red values are not large enough

arrow([0, 0], [725, 350], 0.1, 0.4, 0.1, ...)

There exists a plots:-arrow and a plottools:-arrow (the one you use here given the loading order of the packages).
Depending on the situation I often prefer to use the first one for the width of the arrow and the length of its head can be defined in "relative" dimensions.
For instance width=[0.1, relative] means that the with of tha arrow is 10% of its length.

Here is a proposal for your arrow?
I suppressed the axes and the view option to focus on the arrow itself.

with(plots):
l1 := arrow(`<,>`(0, 0), `<,>`(725, 350), width=[0.02, relative], head_length=[0.1, relative], head_width=[3, relative], color = red, legend="vector v"):
display(l1, axes = frame, axes=none);

Arrow.mw

Rewrite the matriw in a more abstract form
 

MyKGff := Matrix(5, 5, 
  {
    (1, 1) = A+(1/63368)*B, 
    (1, 2) = -(1/63368)*B, 
    (1, 3) = -(1/39605)*B, 
    (1, 4) = 0, 
    (1, 5) = 0, 
    (2, 1) = -(1/63368)*B, 
    (2, 2) = (1/63368)*B+C, 
    (2, 3) = (1/39605)*B, 
    (2, 4) = -C, 
    (2, 5) = 0, 
    (3, 1) = -(1/39605)*B, 
    (3, 2) = (1/39605)*B, 
    (3, 3) = (8/198025)*B, 
    (3, 4) = 0, 
    (3, 5) = 0, 
    (4, 1) = 0, 
    (4, 2) = -C, 
    (4, 3) = 0, 
    (4, 4) = C+(1/63368)*F, 
    (4, 5) = -(1/39605)*F, 
    (5, 1) = 0, 
    (5, 2) = 0, 
    (5, 3) = 0, 
    (5, 4) = -(1/39605)*F, 
    (5, 5) = (8/198025)*F
  }
);

Then (Maple 2015.2)

PseudoInverse := CodeTools:-Usage( LinearAlgebra:-MatrixInverse(MyKGff, method=pseudo) ):
memory used=28.38MiB, alloc change=68.00MiB, cpu time=312.00ms, real time=314.00ms, gc time=24.90ms

To recover the pseudo inverse with original data do

Rewrite := [
             A = (D__pile*W[1] . E__c)/3000, 
             B = (D__pile*W[2] . E__c)*89^(1/2), 
             C = (D__pile*W[3] . E__c)/1000, 
             F = (D__pile*W[4] . E__c)*89^(1/2)
           ]:

eval(PseudoInverse, Rewrite):

InverseMatrix_mmcdara.mw

With Maple 2015.2

restart:
A__1 := A__0 + k(e)*a__1;
diff(phi(A__1), a__1);
                        A__0 + k(e) a__1
                 D(phi)(A__0 + k(e) a__1) k(e)
convert(%, diff);  

which reads "derivative of phi with respect to its argument, evaluated at point A0+k(e)*a1, times k(e)"

This edited version of my first answer (which I removed) owes a lot to @tomleslie  whose answer allowed me to better understand what you wanted.

I suggest you first remove all the spaces from your file.
These commandas hafe to run from a terminal window

sed -i bk  's/ //g' OriginalFile.txt   # Mac OSX
sed -i 's/ //g' OriginalFile.txt       # Windows & LINUX

(on my Mac [OSX] the Maple command 

F   := cat(currentdir(), "/OriginalFile.txt"):
cmd := cat("sed -i bk  's/ //g' ", F):
ssystem(cmd):

returns an error, which is the reason I ran it from a terminal windox)

Now the original file (load it here  OriginalFile.txt) now doesn't contain any space character (look to the meaning of the "-i" option on the web).
This Maple code represents the content of your file as a table with as many entries as data blocks
WithSpacesRemoved.mw

The expression

k/(1-exp(-Pi*k))

cannot be the Fourier transform of FTI for it doesn't contain epsilon.
So I guess this expression could be the Fourier transform of L(0, k):

L := (k, epsilon) ->  -(int(exp(-I*k*(x+I*epsilon))/sinh(x+I*epsilon)^2, x = -infinity .. infinity))/(2*Pi):

, which is not:

L(k, 0)
                             infinity
                           - --------
                                Pi   

Numerically:

M := (k, epsilon) ->  evalf(Int(exp(-I*k*(x+I*epsilon))/sinh(x+I*epsilon)^2, x = -infinity .. infinity, method=_d01amc))/(2*Pi):
M(1, 0)
Error, (in evalf/int) NE_QUAD_NO_CONV:
  The integral is probably divergent or slowly convergent.

matrixplot(x(t), t = 0 .. 4) can't work.
Either you do this matrixplot(x(1)) for a specific value t=1 or, which I guess is closer from what you want:
 

restart
with(plots):
with(LinearAlgebra):
x := t -> `<|>`(`<,>`(.9582+1.9874*t-2.059*t^2+2.6184*t^3-3.7508*t^4, 1.11935+1.4898*t-.6319*t^2+1.037*t^3-1.4695*t^4)):

animate(
  matrixplot, 
  [
    x(t), 
    heights=histogram,
    axis[1]=[tickmarks=[1.5=1, 2.5=2]],
    axis[2]=[tickmarks=[1.5=1]]
  ], 
  t=-1..2
)

matrixplot.mw

As your question is not clear, here are other possibilities OtherWays.mw

I think that Maple (I use Maple 2015) cannot find a formal solution.
Nevertheless one can do a few things:

• simplify the solution by the change of variable xi=2*h/a,
• account for the symmetry wrt xi to reduce the search to [0..+infinity),
• observe that there is a solution iif abs(L) >= 2*abs(a),
• observe that L=2*a ==> xi = 0 ==> h = 0,
• account for a scaling of the solution:
  If h(a, L) is the solution for some couple (a, L) than h(s*a, s*L) = s*h(a, L) for any real s,
• finallly use fsolve instead of solve.

See here
A_few_ideas.mw

A squared matrix which contains only one "1" on each line and each column and "0' everywhere else, is named a LHD matrix (LHD = Latin Hypercube Design).
I will say "LHD of 1" for reasons that will appear later.
Let N the dimension of such a matrix.

There are obviously N!  "LHD of 1", thus 6 "LHD of 1" in your case.

Now I'm going to replace, for one arbitrary "LHD of 1" among these six, some "0" by "1"

1. Change nothing... 1 possibility (meaning "keep the original LHD")
    
2. Build a single row with 2 "1" by starting f from the  "LHD of 1"
   This row can be chosen in 3 different ways and the "1" can be put in 2 different places.
   This lead to 3*2 = 6 different possibilities for a given "LHD of 1".
    
3. Build two rows 2, both with 2 "1" by starting from the  "LHD of 1".
   These two rows can be chosen in 3 different ways (rows (1, 2), (1, 3) or (2, 3)).
   On the first row one can put the second "1" at 2 different places. Let c the number of the column this "1" is placed.
   There are two cases to consider:
   1. The second chosen row contains a "1" in column c: then there are 2 places left to place a second "1" on this row.
   2. The second row contains a "0" on column c: if one puts a "1" on column c, this column will necessarily contains 3 "1".
      So there is only 1 place left on the second row where a "1" can be placed.
   3. Finally there are 3*(2+1) = 9 possibilities to build two rows with 2 "1" per  "LHD of 1".
       
4. Finally, build three rows with 2 "1" by starting from the  "LHD of 1".
   If one observe that this gicves a "LHD of 0", there are only 1 way to do that for a given "LHD of 1".
    

If we sum up these results, there are  6 * (1 + 6 + 9 + 1) = 102  matrices that meet your requirements

As they are probably different interpretations of your question, here is one
An_Idea.mw

The trick I used here is the transformation of each "Rieman" polygon into a prism by using plottools:-extrude.

Another method could be to construct  the piecewise function which represents the envelope of these "Rieman" polygons and make it turn (of course the rotating figure would then be a flat figure).

I think that exact integration is possible only if y = 0

Look at this

restart:
J := (y, tau) -> -2/3 + int(sqrt(x)/(exp((x-y)/tau ) + 1), x=0..infinity):
J(0, t) assuming t > 0;
J(0, t) assuming t < 0;
J(1, t) assuming t > 0;

In case you would be interested in numerical values only, you can do this

J := (y, tau) -> evalf(-2/3 + Int(sqrt(x)/(exp((x-y)/tau ) + 1), x=0..infinity, method=_d01amc)):
J(1, 1)
                          0.7297086143