i am almost sure that those answering (if having time and interest) do not even need points or awards but mostly may prefer something like "thank you, that helped me", if it is honest

  restart; Digits:=20;
  k/Pi^4; subs(k=104,%); a:=evalf(%); b:=1/a;

  identify(a, extension=[seq(Pi^(j),j=-4..4)]);

                                 104
                                 ---
                                   4
                                 Pi

  identify(b, extension=[seq(Pi^(j),j=-4..4)]);

                                   4
                                 Pi
                                 ---
                                 104

Convert to listlist or array (the last suggested through the help pages) works for me.

But it forgets the floating point nature: 2.0 becomes an integer 2

Have not stepped through the problem in detail, but you do not need
to approximate the error function on your own (it makes the function
non-analytic) and there is no need to use Simpson (playing with your
approach I see a 'discontinuity' at ~ 1.6 for your task, but do not
want to look closer).

After converting to rationals you want the following:

Int(169858969913624935089195720356683/26813033068808184917956157510272*
  (1-erf(5/12*(12/5-u)/u))/u*exp(-1/9*(u-1/2)^2+25/144*(12/5-u)^2/u^2),
     u = 125/942 .. 175/314);

Using simplify(%); expand(%); combine(%,exp); normal(%); evalf(%);
gives the value in the title (astonishingly not 64 digits ...)

hm ...

You said you want the images ...

... but if you want the complete 'sheet' why do you want to run
it from outside and not some sheet itself?

Note that you do not have a sheet (to be saved), if you use only
the kernel itself ... but you want the userinterface to be shown
from it (I am not aware that one would do so with other systems,
only SAP had something similar called 'batch input').

It would be far more simple to enter the commands directly than
to use the cryptic syntax for OpenMaple ...

May be you can fake some protocol, save that and reload it.

Or do you want to 'hack' the 3rd party solution? There should be
a reason, why the software allows you no deeper access ... Do
you have the source code?

If not you may also contact the software vendor, whether (s)he
offers a way to achieve your needs.

what should that be good for, why you want it that way?

A better method is to use RootOf (and to multiply by 1+k^2), but it gives no 'answer'
even after applying 'allvalues'.

You may want to interprete it as: "the system can not provide a satisfactory answer
automatically in terms of its known functions", so it is quite likely that there is none
(even if there may be a solution [depending on conditions for k]).

Of course t=0 will do, but for others ...

 exp(-t)-cos(t)+sin(t); # k=1
 fsolve(%,t,complex,avoid={t=0.}); 

                 -6.5025973236926 - 6.8491704362806 I

I think the question is a numerical variant for the the other thread
http://www.mapleprimes.com/forum/solve0 by the same poster and there
is no reason (for me) why a 'symbolic solution exists' or why any CAS
should find one.

For the problem here: he is working only with 6 significant Digits,
and must be aware for the consequences (besides possible typos ...).

It is not even clear whether the small constants are to be considered
as zero, a variant to acer's post is:

  Digits:=14;
  MmaTranslator:-FromMma(p); 
  fnormal(%);                       # clean up
  lhs(%)-rhs(%): convert(%,rational);
  subs( t= 200000/76603*s, %);      # to have cos(1*s) and sin(1*s)
  q:=%;

The above cleaned up the input and normed it a bit, we want zeros for
q depending on s and for s=0 we have one (without cleaning up through
fnormal it comes to acer's result, the taks is symmetric w.r.t. the
'y-axes'). Due to periodics another solution is in s + a, a=76603*Pi
for the above presentation.

It may be better to give some information about the actual task which
rised the question(s).

Why there should be a symbolic solution? I see no reason for that.
But it is not difficult to find the numerical ones:

Calling R the initial expression you can combine tau and h into t,
RR:=subs(tau=1,h=1,R), M:=subs(t=I*t,RR) following Alex' way, now

  map(Re,M);
  plot(%, t= - 2.3..2.3):
  plot(1, t= - 2.5 .. 2.5, color=grey):
  plots[display](%%,%);

shows the roots for M=1, they are symmetric w.r.t. vertical axis:

  'eval(M,t=-t)' - 'M': '%'= simplify(%);

                           M|       - M = 0
                            |t = -t

  'evalc( Im(M))': '%'=%;
                           evalc(Im(M)) = 0

and

  fsolve(M=1,t=0.5, fulldigits); eval(RR, t= %/I); evalf(%);
  fsolve(M=1,t=2.1, fulldigits); eval(RR, t= %/I); evalf(%);

finds and confirm the roots for RR - 1, they are given through
0.56880852623983 and 2.2446321162676 (Maple12, 14 Digits).

Knowing nothing I googled and ... besides that I got the impression that the ODE is not always solvable in closed form was led to Zwillinger "Handbook of Differential Equations (1997)", who in 44.1.1, p.180 refers to Murphy, G. M. "Ordinary Differential Equations and Their Solution" (1960) at pp. 23 and 25. May be time for OP to visit a library ...

in both Maple 11 and 12 I get exp(mu2(r)) by just pasting the expression, not even a simplify is needed (classical interface)

best would be you upload your sheet

I do it the following way having MS Word (and having always activated non-printable signs):

- copy and paste into a Word document (may be with font=courier for better reading)
- then search for  "^p> "   and replace it by "^p" (all occurencies)
- and then copy and paste to Maple (I only use classical interface)

The above means to replace a paragraph sign ^p followed by a > and a blank is
replaced by the paragraph sign only (which is the hard line break)

Literally taken for me the question also makes no sense and as an
ordered finite set is nothing but the according interval of the
Naturals for me that would sufficient (and use substitutions only).

A different source (as I do not have Knuth) and as C program is in
the beautiful online book Jörg Arndt, "Algorithms for programmers",
Chapter 10 at www.jjj.de/fxt/

Have not checked, whether it is Algorithm L.

Thx for the info, but I do not know anything about physical Chemistry :-)
Still I would say, that *your data can not fit the model*, same reasoning.

However I modified your model using 2 location parameters and 

  G := x -> - c
  -a^2*exp((x-mu)*r^2)/(1+exp((x-mu)*r^2))+b^2*exp(x*r^2)/(1+exp(x*r^2))^2;

                        2               2        2        2
                       a  exp((x - mu) r )      b  exp(x r )
        G := x -> -c - -------------------- + ----------------
                                         2                2  2
                       1 + exp((x - mu) r )   (1 + exp(x r ))

From your data and the limiting behaviour at +-infinity of the original
model it follows, that L must be positive and A is negative:

  '[limit(EG(x), x=-infinity),limit(EG(x), x=infinity)]': %: % assuming L<0: 
    '%%%'=%, L<0;
  '[limit(EG(x), x=-infinity),limit(EG(x), x=infinity)]': %: % assuming 0<L: 
    '%%%'=%, 0 < L;

   [      lim         EG(x),      lim       EG(x)] = [A, 0], L < 0
    x -> (-infinity)         x -> infinity


   [      lim         EG(x),      lim       EG(x)] = [0, A], 0 < L
    x -> (-infinity)         x -> infinity

where EG is your original model function, EG:= x -> ((A*exp((2*Pi*x)/L))/
(1+exp((2*Pi*x)/L)))+((B*exp((2*Pi*x)/L))/(1+exp((2*Pi*x)/L))^2);

Trying to fit that always shows a shift off zero: the maximum can not be
achieved close to it, because of A being not 0, since it is responsible
for the tail as well - at that is below -0.02 and certainly lower. 

So I introduced a location parameter mu: some shift along the x-axes for
the "A-term".

Additionally your data at the left side of zero stop at the x-axes.
But I guess that's only true for the provided data and you simply do
not have more. So it seems natural for me to assume data may become
negative as well there. But the limiting behaviour wants them to be
strictly positive.

So I use a level shift parameter c, shifting vertically.

And I came up with the definition for G (close as possible to your model
and the given data):

  EG(x); subs(A=-a^2, B=b^2, L=Pi*2/r^2, %);
  subs(x=x-mu*1,op(1,%)) + subs(x=x-nu*0,op(2,%));
  -c + %;
  G:=unapply(%,x);

Then I get a reasonable fit:

  residuals := map(p -> (G(p[1])-p[2]), data):
  R:= LSSolve(residuals);
  params:=R[2];
  Gfit := unapply(eval(G(x), params),x);
  plot(gfit(x), x=-6 .. 6);
  display([pointplot(data), plot(Gfit(x), x=data[1,1]..data[-1,1])]);

Certainly not superb, but reasonable.

will try tonight ... if data go negative to the right that should bear info on signs of A and L, see above