Changing variable helps:
f := (t-1)^(1/3);
p:=2;
b := 2/p*Int(f*sin(2*Pi*n*t/p), t = 0 .. p);
b1:=IntegrationTools:-Change(b,t=1+u);
bv:=value(b1);
#Check:
evalf(eval(bv,n=3));
evalf(eval(b,n=3));
###
As Carl points out, if you want the real third root you need surd instead.
subs(u^(1/3)=surd(u,3),b1); #replacing with surd
value(%); #Doesn't work, however
#So we do it by using that u^(1/3)=-(-u)^(1/3) when u<0 and the real root is desired
IntegrationTools:-Split(b1,0);
bneg,bpos:=op(%);
bposv:=value(bpos);
bneg2:=subs(u^(1/3)=-(-u)^(1/3),bneg); #Instead of surd
IntegrationTools:-Change(bneg2,u=-t);
bnegv:=value(%);
evalf(eval(bposv+bnegv,n=3));
evalf(eval(bpos+bneg2,n=3));

If you right click on your 2-D Math you have a choice of converting it to 2-D Math or 2-D Math Input (or 1-D Math Input). That makes a difference.

Assuming that you are satisfied with pdsolve doing all the work then do:
pde:=diff(u(x,t),t,t)=c^2*diff(u(x,t),x,x);
pdsolve({pde,u(x,0)=f(x),D[2](u)(x,0)=g(x)});
res:=IntegrationTools:-Combine(%);
eval(res,{c=4,f=(x->1/2*exp(-x^2)),g=(x->exp(-x^2))});
res1:=value(%);
plots:-animate(plot,[rhs(res1),x=-42..42],t=0..10);

Since eq1=0 implies f'''(y)-f'(y)=constant your P2[x] is constant for each x. To get the constant you can use y=0 since it lies in the interval h2..h1 for all x.
Thus change the definition of Vls to
Vls:=Vector([seq(P2[x](0),x=0..1,0.1)]);

Agreeing totally with Carl I tried numerically.
The idea is to solve the 8 equations using dsolve/DAE.
First we need to find initial conditions for the functions in the algebraic equations. This can be done like this, where I'm using your sys and dconds:
eqs0:=subs(t=0,dconds,sys);
indets(%,function);
res0:=solve(eqs0,%);
res01:=remove(has,[res0],I); #Removing imaginary solutions
res:=op(remove(hastype,res01,negative)); #Removing negative solutions
#Alternatively in one line:
res:=op(remove(hastype,[res0],{imaginary,negative}));
#Then proceed with defining dsys as you did and finally
dsol:=dsolve({dsys,dconds} union sys union res,numeric);
plots:-odeplot(dsol,[[t,Pct_Al_i(t)]],0..1000); #Example plot

Notice that
seq(laplace(Dirac(i,t),t,s),i=0..4);
returns
1, s, s^2, s^3, s^4
Notice also that your expression Y[1] (after you redefine it) is a rational function of s:
type(Y[1],ratpoly(anything,s));
and furthermore the degrees of the numerator and denominator are 13 and 9, respectively:
degree(numer(Y[1]),s);
degree(denom(Y[1]),s);
Thus you should expect the appearance of the delta function and its derivatives up to order 4 (=13-9).

Incidentally, what is the rationale for setting I=0 later on?

I don't see any way of finding those times without solving the system of odes. Since you say that your system contains lots of terms it seems unlikely that you can find closed form solutions, thus you need to solve numerically.
Here is a rather trivial example which can be solved symbolically, but I also solve it numerically:

sys:=diff(a(t),t,t)=-a(t),diff(b(t),t,t)=-4*b(t);
dsolve({sys,a(0)=1,D(a)(0)=0,b(0)=0,D(b)(0)=1});
`=`(op(rhs~(%)));
solve(%,t,allsolutions);
#Now solving numerically. I'm using events. The event triggers when a(t)=b(t) and the action is to give c(t) that value until the event is triggered again:
res:=dsolve({sys,a(0)=1,D(a)(0)=0,b(0)=0,D(b)(0)=1,c(0)=0},numeric,discrete_variables=[c(t)],events=[[a(t)=b(t),c(t)=t]]);
res(2); #Example
plots:-odeplot(res,[[t,a(t)-b(t)],[t,c(t)]],t=0..20);
#If you would like the set of times, you could modify to:
res:=dsolve({sys,a(0)=1,D(a)(0)=0,b(0)=0,D(b)(0)=1,c(0)=0},numeric,discrete_variables=[c(t)::float],events=[[a(t)=b(t),c(t)=t]],output=listprocedure);
C:=subs(res,c(t));
{seq(C(i*.1),i=1..200)}; #t=0 doesn't count. It is just the initial value given to c.
#Compare with the exact result:
evalf(seq(Pi/2*(2*i-1),i=1..6));

You should correct vector to Vector in the definition of upsilon.
Your matrix Lambda can be defined like this:
Lambda:= Matrix(2,(i,j)->upsilon[j]*exp(-lambda[i, j]/(Rgas*T))/upsilon[i]);
Notice that your initial assignments to i and j are irrelevant here (and you don't need them). Here they are just used as formal parameters in a procedure.

You can try
dsolve(Eq1);
The result doesn't make me think that you have any chance of solving for the conditions in bcs: Notice the integrals.
But you can solve numerically:
params:={k=1,B=1,h=1};
res:=dsolve(eval({Eq1,bcs},params),numeric,method=bvp[midrich]);
plots:-odeplot(res,[r,f(r)],-1..1);

You can use map without any assumptions since it just does what you ask it to do:
restart;
c:=a<b:
map(x->x/d,c);

I shall ignore phi since it doesn't play any role. Also for convenience I use w, s1, and s2 instead of the indexed versions.
This is a first order linear pde.
restart;
sys:=-(1/2)*(diff(f(phi, s[1], s[2], w[1]), s[1]))*(s[1]-s[2])*s[1]/(s[2]*w[1])+(1/2)*(diff(f(phi, s[1], s[2], w[1]), s[2]))*(s[1]-s[2])*s[2]/((-1+w[1])*s[1])+diff(f(phi, s[1], s[2], w[1]), w[1]) = 0;
pde:=subs(phi=NULL,s[1]=s1,s[2]=s2,w[1]=w,sys);
#Finding the characteristic odes using w as independent variable:
diff(f(S1(w),S2(w),S3(w)),w);
pdeD:=convert(pde,D);
SYS:=seq(diff(cat(S,i)(w),w)=coeff(lhs(pdeD),D[i](f)(s1,s2,w)),i=1..3);
SYS2:=subs(s1=S1(w),s2=S2(w),S3(w)=w,{SYS[1..2]}); #The characteristic odes
res:=dsolve(SYS2);
res2:=simplify(op(eval(res[2],res[1])));
res1:=op(res[1]);
#Now f(S1(w),S2(w),S3(w)) is constant along each characteristic curve, thus
eval(f(S1(w),S2(w),S3(w)),{res1,res2,S3(w)=w});
#is constant in w, so equal to
eval(%,w=1);
#i.e. an arbitrary function g of _C2. Now find _C2 from res1 and res2:
solve({res1,res2},{_C1,_C2});
#Thus solutions are
sol:=f(s1,s2,w)=subs(%,S2(w)=s2,S1(w)=s1,g(_C2));
pdetest(sol,pde);

The problem is that solve cannot find explicit representations for the roots of the polynomial of degree 6 you got in the second case.

One solution is to use fsolve instead of solve, e.g. like this:
p:=2: q:=3:
ek:=unapply(eq,K);
U:=K->[fsolve(ek(K),u,complex)];
U(.5);
animate(complexplot,['U(K)',style = point,symbolsize=10,color="red"],frames=50,K=-1..1,background=unitcircle, scaling=constrained, trace=20 );

An alternative, which is shorter is:
animate(complexplot,[[allvalues(solve(eq, u))],style = point,symbolsize=10,color="red"],frames=50,K=-1..1,background=unitcircle, scaling=constrained, trace=20);

No way. a*b and b*a are kept in the same place in memory.
Consider these two sessions (where I also use collect as suggested by nm):

restart;
g:=a*x+a*y;
factor(g);
collect(g,a);
restart;
g:=a*x+a*y;
collect(g,a);
factor(g);

You can get the full source code for LinearAlgebra:-Permanent like this:
showstat(LinearAlgebra:-Permanent);
showstat(`tools/ClearRememberTable`);
showstat(LinearAlgebra::`Permanent/pminor`); #Notice ::, not :-
#`Permanent/pminor` is local to LinearAlgebra

One obvious change you can make to your code is to replace the inner loop with add so that the double loop becomes:
for i to m do
     rowsum:=add(B[i,j],j=1..s);
     prod := prod*rowsum;
end do;
#
In fact you can replace the double loop by:
prod:=mul( add(B[i,j],j=1..s),i=1..m);
#
Also while at it, it should be an advantage not to call LinearAlgebra:-SubMatrix, but use a more direct approach:
permanentRyser3 := proc (M::Matrix)
      local S, T, B, m, n, s, i, j, sum, prod;
      m, n := op(1, M);
      if m <> n then
               error "expecting a square Matrix, got dimensions %1, %2", m, n
      end if;
      sum := 0;
      S := combinat:-subsets([seq(1 .. m)]);
      S[nextvalue]();
      while not S[finished] do
             T := S[nextvalue]();
             s := numelems(T);
             B:=M[1..m,T];
             prod:=(-1)^s*mul( add(B[i,j],j=1..s),i=1..m);
             sum := sum+prod;
      end do;
      expand((-1)^m*sum)
end proc:

Notice that your map version first evaluates a seq, then uses map. Double work. So a priori the second version ought to be better.