I do not understand what do you mean by the sentence:

I want to solve the system basically considering p,q and their derivatives as constants.

If you say, e.g., that diff(p(x,t),x$3)=k1 in your system, then p(x,t) is not generically a constant:

pdsolve(diff(p(x,t),x$3)=k1);
                       3
                   k1 x                2
         p(x, t) = ----- + 1/2 _F1(t) x  + _F2(t) x + _F3(t)
                     6

Is there a typo above? As a function of s, (1+s*D)/(A*s^2 + B*s + C) = (1+s*D)/(s^2 + 2*damp*wn*s + wn^2) only if A=1.

Your actual input contained things like a[0, t]? May be that you intend to use diff_table? (see ?PDEtools/diff_table)

The general rule is that active commands like sum are for explicit computations, like perform the sum when possible. While inert commands like Sum are more suitable for algebraic manipulations, like in the calculation you intend to reproduce. 

Usually it is not very easy. On one hand, computers do not work like human minds. On the other hand, program designs reflect the mindset of their developers, that can be quite different of the users' mindset.

You can get lots of information if you wish. Uncomment one or more of the commands, and play with the numeric values:

#trace(dsolve):
#infolevel[dsolve]:=2:
#printlevel:=10:
dsolve(diff(y(x), x)-x*y(x)+exp(-y(x)));

As defined, nIntConv produces nested procedures and it makes sense only for a natural (positive integer) argument, i.e. 1, 2, etc, not for a symbolic argument n:

rect := t -> piecewise(t < -1/2,0,1/2 < t,0,1):

intConv := (f1, f2) -> t -> int(f1(t-tau)*f2(tau),
tau = -infinity .. infinity):

nIntConv := proc (n) if n = 0 then RETURN(intConv(rect, rect)) 
end if; RETURN(intConv(nIntConv(n-1), rect)) end proc:

nIntConv(n);
Error, (in nIntConv) too many levels of recursion

So, this error message is not an issue of limit.

It can be seen more simply:

w*a+w[j];
                              w a + w[j]

algsubs(w=beta*w[j],%);

Error, (in collect/series) invalid expression for series

And can be traced here:

series(w[j], w);

Error, invalid expression for series

Most likely is a mess involving indexed names. So, I would avoid them and use something else like:

w*a+w_j;
algsubs(w=beta*w_j,%);
                              w a + w_j

                           a beta w_j + w_j

You can do something like:

with(DEtools):
sys := {diff(x(t),t) = v(t), diff(v(t),t) = -.4*v(t)-1/2*sin(x(t))}:

DEplot(sys,[x(t),v(t)],t=0..20,[[x(0)=1,v(0)=0],[x(0)=0,v(0)=-1],
[x(0)=1+Pi,v(0)=0],[x(0)=Pi,v(0)=-1]],arrows=none);

And add more curves by adding more ICs.

You have just done it at pasting them in your post. Clearly,  x[2] (name[something]) is an indexed name, while `#msub(mi("z"),mn("1"))`  (`#mathmlcommand(something);`) is a "literal" name.

You can do e.g :

`#msub(mi("z"),mn("1"));`:
lprint(%);

`#msub(mi("z"),mn("1"));`

On one hand you pose a 2nd order ODE for t(x), and on the other hand your BC1 already gives t(x)...

In any case, c=1 makes  c + (c-1)*erf(x) = 1 for all x.

This may be closer:

writedata[append]('terminal',[-0.123456-0.123456*I, 
0.1334423423*10^(-15),3], float, proc(f,x) fprintf(f,`%.5Zf`,x) 
end proc);

-0.12346-0.12346I
1.334423423e-16
3

I think that in this case, it is more an issue of flags that control the route of the simplification. See, here:

showstat(`simplify/do`,1..2);
`simplify/do` := proc(s)
local f, i, inds, r, rnormal, symb_mode, new_inds, kernels, 
new_kernels, new_r, new_simp, op0fun, tried_simplify_constant;
   1   if 1 < nargs and member(args[nargs],`simplify/options`) then
   2     symb_mode := args[nargs]
       else
         ...
       end if;
       ...
end proc


`simplify/options`;
       {ask, encode, proviso, symbolic, remove_multiplicities}

The evaluation of sin(0) as 0 occurs because this information is contained in the permantent cache of sin. In a fresh session: 

op(4,eval(sin)); 
                                    Pi
  Cache(512, 'permanent' = [0 = 0, ---- = 1, I = sinh(1) I,
                                    2

                                       1/2
                                Pi    2
        -infinity = undefined, ---- = ----, infinity = undefined,
                                4      2

                1/2
         Pi    3              Pi
        ---- = ----, Pi = 0, ---- = 1/2])
         3      2             6

This normalization is done, I understand, for efficiency reasons in computation. For algebraic manipulations, where immediate computation is not required, it may be more convenient to use the inert form of the function call:

%sin(0);
                               %sin(0)

value(%);
                                  0

For instance exp(x)+3/2*x^2:

x:=c[2]*e/c[1]+(2*c[3]*c[1]-3*c[2]^2)*e^2/c[1]^2+(3*c[4]*c[1]^2-
10*c[2]*c[3]*c[1]+8*c[2]^3)*e^3/c[1]^3:

simplify~(series(exp(x)+3/2*x^2,e=0,3));
                                            2
                 c[2]     2 c[3] c[1] - c[2]   2      3
             1 + ---- e + ------------------- e  + O(e )
                 c[1]                2
                                 c[1]

There has been several threads here about indenting code and use of external editors. See e.g. this one.