I was asking what dsys was assigned, you left that out of the steps shown. Probably it is

dsys := { ec1, ec2, <initial conditions> };

where <initial conditions> are just that, maybe alpha(0) = 0, beta(0) = 0. 

A minor problem is that tau[1] and tau[2] are not assigned numeric values.  They could be declared as parameters to dsolve, but at some point they would have to get numeric values before integrating. The more serious issue is that dsolve cannot convert your set of differential equations to a first order system.

The usual way to combine plots is with ?plots[display].  For example,

y1 := x+4:
y2 := 2*x+1:
p1 := plot(y1, x=0..5, color=red):
p2 := plot(y2, x=0..10, color=blue):
plots:-display(p1,p2);

Don't worry about the code appearing to be "cut off"; we can copy the entire block without a problem.

A nice way to solve the problem is to extend your system so that the integration is performed by dsolve.  To do that, just add the differential equation diff(tot(t),t) = S(t) + R(t), with tot(0) = 0.  Thus,

with(plots):

deqs := {diff(S(t), t) = t^2
        , diff(R(t), t) = t
        , diff(tot(t),t) = S(t)+R(t)
       }:

ini := {S(0) = 1, R(0) = 1, tot(0)=0}:

ode1 := dsolve(deqs union ini, numeric
               , 'events' = [[t = 600, S(t) = 10^7], [t = 800, S(t) = 5*10^7]]):

odeplot(ode1, [t, S(t)+R(t), color = blue, thickness = 2], 0 .. 1000, numpoints = 1000); 
ode1(1000);
[t = 1000., R(t) = 500001.000000000, S(t) = 0.212666666666667 10 ,

                                 11
    tot(t) = 0.469666682666667 10  ]

What do you plan to do with them?  You might consider writing a simple module that provides some conversions.  For example,

krylov := module()
option package;

global K1,K2,K3,K4;

export To, From, Simplify, Evalf;

local ModuleLoad, removeImaginary;

    ModuleLoad := proc()
    global `convert/tokrylov`, `convert/fromkrylov`, `simplify/krylov`
        , `evalf/K1`, `evalf/K2`, `evalf/K3`, `evalf/K4`;

        `convert/tokrylov`   := To;
        `convert/fromkrylov` := From;
        `simplify/krylov`    := Simplify;
        `evalf/K1`           := x -> Evalf(K1(x));
        `evalf/K2`           := x -> Evalf(K2(x));
        `evalf/K3`           := x -> Evalf(K3(x));
        `evalf/K4`           := x -> Evalf(K3(x));

    end proc;

    To := proc(ex)
        eval(ex, [cos = K1-K3,
                  sin = K2-K4,
                  cosh = K1+K3,
                  sinh = K2+K4
                 ]);
    end proc:

    From := proc(ex)
        eval(ex, [K1 = 1/2*cosh + 1/2*cos,
                  K2 = 1/2*sinh + 1/2*sin,
                  K3 = 1/2*cosh - 1/2*cos,
                  K4 = 1/2*sinh - 1/2*sin
                 ]);
    end proc:

    Evalf := proc(ex)
        evalf(From(ex), _rest);
    end proc;

    Simplify := proc(ex)
        # a very basic simplification, Kx(I*z) --> Kx(z)
        evalindets(ex
                   , 'specfunc(anything,{K1,K2,K3,K4})'
                   , fx -> op(0,fx)(removeImaginary(op(1,ex)))
                  );
    end proc:

    removeImaginary := proc(ex)
        if ex :: imaginary
        or ex :: `*` and ormap(type,ex,imaginary)
        then
            return ex/I;
        else
            return ex;
        end if;
    end proc:

    ModuleLoad();

end module:

with(krylov):

Simplify(K1(I*x));
                                     K1(x)

convert(exp(2*x),'compose,trigh,tokrylov');
                     K1(2 x) + K3(2 x) + K2(2 x) + K4(2 x)

evalf(K1(Pi/2 + K3(3)),30);
                        303.220352664355942804059979113

This is just something I threw together in a few minutes; you'd want to rewrite it to do something useful

You need to reformulate the equation for a(t) so that it is not in terms of a(t).  The dsolve command would normally do that, but the piecewise hinders it. You also need to replace time with t. The following should do what you want

m := 5; F0 := 10; a0 := 5;
deqs := { m*a(t) = piecewise(F0/m*t < a0, F0*t, 0)
          , v(t) = diff(x(t), t)
          , a(t) = diff(v(t), t)
        }:
ics := {x(0) = 0, v(0) = 0};
integ := dsolve( deqs union ics, numeric):

plots:-odeplot(integ, [t,x(t)], 0..3);

Note that for the units to be consistent, F0 has dimensions of force/time.

Followup

I see that does not work past t=2.5.  One way around that is to use dsolve event handling. While this can be a bit complicated, in this case we can handle it simply by using a discrete variable (delta) whose value changes from 1 to 0 when F0/m*t-a crosses zero. A better model might pay attention to the direction---not an issue here with the given equations and initial conditions.

deqs := { m*a(t) = F0*t*delta(t)
          , v(t) = diff(x(t), t)
          , a(t) = diff(v(t), t)
        }:

ics := {x(0) = 0, v(0) = 0, delta(0)=1};

# Use a single event that is triggered when the acceleration
# reaches the threshold.  At that time, change delta to 0
# which causes the effective acceleration to go to 0.
evts := [[F0/m*t-a0=0, delta(t)=0]]:

integ := dsolve( deqs union ics, 'numeric'
                 , 'events' = evts
                 , 'discrete_variables' = [delta(t)::float]
               ):

plots:-odeplot(integ, [t,v(t)], 0..3);

The error message indicates what needs to be done.  In the assignment to P, change each x to x(t).  Alas, dsolve won't give a symbolic solution, but you can use dsolve,numeric.

integ := dsolve({eq,ics}, numeric):
integ(1);
plots:-odeplot(integ, [t,x(t)], 0..1);

@Chef-Chaudard It would be helpful if you showed both what you did and what Maple returned.  Are you sure you replaced the single forward quotes ('...') with double quotes ("...")? The double quote is a single symbol, not two single quotes.

The format operand to fprintf should be a string; a name may work but only if Maple can parse it, and here it won't parse the $, which is an operator. So enclose the format operand in double quotes.

Are you asking whether the indexing operation can be overloaded?  The answer is yes, though it may not be the best way to go.  Here is an example, I've modified your MyType definiition to be an inert function rather than a list. It could be done with a list, but I'm not seeing a good reason to do so:

restart;
MyPkg := module()
option package;

export `?[]`;
local ModuleLoad;

    ModuleLoad := proc()
        TypeTools:-AddType('MyType', 'specfunc(anything,MyType)')
    end proc:

    `?[]` := proc(base,indx)
    local i;
        if base :: MyType then
            op(indx,base);
        else
            :-`?[]`(base,indx);
        end if;
    end proc;

    ModuleLoad();

end module:

with(MyPkg):
T := MyType(a,b,c);
T[1];
                            a

It depends on what you mean by the source code. If you create and distribute a Maple archive (.mla file), the user can run the procedures and doesn't have access to the actual source code, which includes comments and preprocessor macros.  However, the user can always view the interpreted Maple code.

Maple immediately simplifies fractions, so there is no nice way to do what you want.  As a workaround, you could enclose the numerator and denominator in the empty function (``) and then later call expand. For example,

(**) F := ``(4)/``(6);
                                     (4)
                               F := ----
                                     (6)

(**) expand(F);
                                  2/3

Context would be helpful. What command are you using?

You are correct, the attribute is not copied.  It isn't clear whether that is a bug or by design.   You can access the attribute table by doing

op(op(5, G));
table([0 = ["graphAtt" = 1],
                                          -16
    1 = ["draw-pos-fixed" = [0.61232340 10   , 1.], "vertexAtt" = 1],
    2 = ["draw-pos-fixed" = [0.95105652, 0.30901699]],
    3 = ["draw-pos-fixed" = [0.58778525, -0.80901699]],
    4 = ["draw-pos-fixed" = [-0.58778525, -0.80901699]],
    5 = ["draw-pos-fixed" = [-0.95105652, 0.30901699]],
    (1, 2) = ["edgeAtt" = 1],
    "draw-default-style" = "draw-pos-fixed"
    ])

You are using a matrix, not a Matrix.  The former is an older data structure.  You can substitute into a matrix, but you have to eval it:

subs(a = 1, eval(A));

You could also do

eval(A, [a=1, b=2]);

However, you probably want to use a Matrix. Also, rather than linalg (which does use the matrix type), you should consider using LinearAlgebra, which uses the Matrix type.

V[1] evaluates to the first entry in V, in this case, xt[1, 1]*xt[1, 2]*alpha. That expression is a product.  When select is applied to a product, the returned expression is (generally) a product. If no item matches, 1 is returned.  For example

select(`=`, a*b, 2);
                                     1

To understand why 1 is returned, consider

(x,y) := selectremove( some_predicate, some_product);

You will then find that some_product = x*y, as expected, regardless the predicate. For a similar reason

mul(x, x=1..0) --> 1