You can also do

evalf[100](sqrt(x))

I agree with the original poster, there should be an inverse to CoefficientList/CoefficientVector in PolynomialTools. Creating one is easy enough, but it should really be in the PolynomialTools package. Otherwise less sophisticated users will do so with sum or even term by term in a do-loop.

Using the solution for hyperbolic motion (p. 323) try the following:

eqhyperb := (x+c^2/a)^2 - c^2*t^2 = c^4/a^2:
sols := solve(eq, x);
sol := sols[1]:

eqphoton := c*(t-t0):

asymp := convert(asympt(sol,t,1),'polynom');

asol := solve(identity(asymp = eqphoton,t),[a,t0]);

#check
plot(subs(asol[], a=3, c=1, [sol,eqphoton]), t=0..3);

Unzip the file and add the name of the directory to libname.  On my linux box I did

$ unzip Gravitation.zip
$ cd Gravitation
$ ls
Documentation.mw  Gravitation.mw  maple.ind  maple.lib
$ pwd
/home/joe/maple/lib/Gravitation
$ maple11
libname := "/home/joe/maple/lib/Gravitation", libname:
with(Gravitation);
[aboutGravitation, addObjects, areEqual, createMetric, createObject,
    createTetrad, deriveChristoffel, deriveEinstein, deriveKretschmann,
    deriveRicci, deriveRicciScalar, deriveRiemann, deriveSpinConnection,
    deriveSpinCurvature, deriveWeyl, multiplyObjects, subtractObjects]

Actually, what we would really do is simply

plot(sin, 1..12);

which is shorter and nicer than what is required in Matlab.  No use plotting vectors/lists if you can use Maple's adapative plotting capability. Just because Matlab is restricted to plotting vectors doesn't mean that that is the best way to go.

Not being particularly adept at evaluating integrals, I took a backwards approach and used Maple's identify function to guess the expression for particular values of a.  It was more amusing than enlightening:

# use the half-integral
J2 := a -> Int(exp(-a*x^2)*sech(x/2)^2, x=0..infinity):

guess := proc(a)
local y,id;
    y := evalf[30](J2(a));
    id := identify(y);
    if not id::float then
        'J2'(a) = id;
    end if;
end proc:

for a in [0,1,sqrt(Pi)] do
    guess(a);
end do;
                                   J2(0) = 2

                                       3/5      3/5
                                   10 7    ln(3)
                           J2(1) = -- -------------
                                   49        3/8
                                        ln(2)

                                        3/5        8/5
                             1/2    10 7    Zeta(5)
                        J2(Pi   ) = -- ---------------
                                    49         8/7
                                          ln(3)

You might try, instead, using animate.   I think the result would be more useful.

with(plots):
animate(plot, [x^n, x=-1..1], n=1..20);

fopen("/usr/local/mydir/examp/example.dat");

read "/home/user/example/mycode.mws";

Actually, that last one is questionable----normally one doesn't use the `read' command to directly read a .mws file, those files usually contain a classic style Maple worksheet.  If that file really contains Maple script code, than you might consider renaming it to mycode.mpl.

See the Functional Operators help page (?->):

f := z -> z + 7+I;

My first thought was to use degree in a call to remove, alas, that won't work because your expression is not a polynomial.  Here is one approach.

> f := series(add(k*x^((k-1)/2),k=0..8),x,3);          
                      1/2            (3/2)      2      (5/2)      3
          f := 1 + 2 x    + 3 x + 4 x      + 5 x  + 6 x      + O(x )

> remove(m -> match(m=k*x^p,x,'s') and subs(s,p)>2, f);
                         1/2            (3/2)      2      3
                  1 + 2 x    + 3 x + 4 x      + 5 x  + O(x )

# The structure is not a series, so we can use this trick to eliminate the big-O term
> eval(%,O=0);                                       
                             1/2            (3/2)      2
                      1 + 2 x    + 3 x + 4 x      + 5 x

A nice way to do this is to extend the procedure so that the returned module includes a 'wave' field that is a procedure that returns the value at a given time.  Here is what I did (I'm not confident I computed the waveform correctly, you can fix that). I took the liberty of fixing a bug in the previous (the X should be declared local to the constructor). 

Note: I fixed a bug in the procedure (I forgot to add the dc term, a0).  I combined the plot of the original data with a plot of the fourier waveform, they compare well.

myFourier := proc(period :: positive
                  , amps :: list
                  , numCoeffs :: posint
                  , $
                 )
local t,X;
    module()
    option record;
    export a0,aAreas,bAreas,aCoeffs,bCoeffs,wave;
        a0 := `+`(amps[])/period;
        (aAreas,bAreas) := seq(Matrix(nops(amps),numCoeffs,
                                      (i,j) -> evalf(amps[i]*X(2*(i-1)*j*Pi/period))
                                     )
                               , X in [cos,sin]);
        (aCoeffs,bCoeffs) := seq(Vector[row](numCoeffs,
                                             (i) -> `+`(convert(X[1..nops(amps),i],list)[])*2/period
                                            ),X in [aAreas,bAreas]);

        wave := unapply(a0 + add(aCoeffs[i]*cos(2*Pi*i*t/period) + bCoeffs[i]*sin(2*Pi*i*t/period)
                                 , i = 1..numCoeffs), t);
    end module
end proc:

T := 10:
ys := [14,18.7,9,4.1,6.7,6,6.3,8.4,4,2.9]:
n := nops(ys):
ts := [seq(i/n*T, i=0..n-1)]:

wave := myFourier(n,ys,5);
wave:-a0;
wave:-aAreas;
wave:-bAreas;
wave:-aCoeffs;
wave:-bCoeffs;
p1 := Statistics:-LineChart(ys,xcoords=ts):
p2 := plot(wave:-wave, 0..10):
plots[display](p1,p2);

Another solution is to use create a list (or Vector) and use Statistics[LineChart].  It assumes the domain is the natural numbers (1,2,3,...) though you can provide it with the actual coordinates using xcoords:

For example

X := [seq(-20..20)]:

Statistics:-LineChart(map(x->x^2, X), xcoords=X);

A nice feature of LineChart is that it includes, by default, the points in the plot.

Here's a completely symbolic approach,

restart;
OP := P-O:
OQ := Q-O:
QP := P-Q:
OR := R-O:
QR := R-Q:
eqs := {NULL
        , a1*i + b1*j = OP
        , a2*i + b2*j = OQ
        , OR = 1/3*OP
       }:
elim := eliminate(eqs, {O,P,Q,R}):

# express QR in terms of i and j
collect(subs(elim[1],QR),[i,j]);

                                  (-a2+1/3*a1)*i+(-b2+1/3*b1)*j

N := 5:
assume(seq(T[i]<T[i+1], i=0..N)):
is(T[0]<T[N+1]);
                         true