If I understand correctly, you want to both animations built simultaneously, rather than sequentially. To do this, remove the option  insequence = true . Also the parameters  b  and  J  should be specified in advance.

Example:

b:=1: J:=2:

tau[1] := h-1+(4/3)*(1/sqrt(3)-h+1)*sqrt(3)*(sqrt(3)*((1/4)*J^2-z^2)/(J^2*b));

tau[2] := (4/3)*sqrt(3)*((1/4)*h^2-z^2)/(h^2*b);

plot2 := plots[animate](plot, [tau[1], z = -(1/2)*J .. (1/2)*J, color = red, legend = shear*stress], h = 1 .. 1+1/sqrt(3)):

plot3 := plots[animate](plot, [tau[2], z = -(1/2)*h .. (1/2)*h, color = blue, legend = shear*stress], h = 0 .. 1):

plots[display]([plot2, plot3]);

Should be

factor(x^2-2*a*x+a^2);

                                                     (a - x)^2

Addition: In standard interface in 2d math you can do a space between  a  and  x . Maple regards a space between two symbols as a multiplication sign.  ax  in your code is just a new symbol, not  a*x . A space between a number (numeric type) and a symbol you do not necessarily have to do. For example in 2d math,  2a  is the same as 2 a  or  2*a 

V_matrix:=<seq(V[i], i=1..100)>;

You can solve the equation numerically for specific values of parameters and an initial condition.

Example:

eq:=diff(y(x),x) - (Q - x*p0*(exp(alpha-beta*y(x)))/(1+exp(alpha-beta*y(x))))^2=0:
sol:=dsolve({eq, y(0)=0}, numeric, parameters=[Q,p0,alpha,beta] );
sol(parameters=[1,2,3,4]);
plots[odeplot](sol, [x,y(x)], x=0..5);

 Addition. `Series` solution (above) fits only near x=0. Compare: 

restart;
eq:=diff(y(x),x) - (Q - x*p0*(exp(alpha-beta*y(x)))/(1+exp(alpha-beta*y(x))))^2=0:
P:=[Q,p0,alpha,beta]:
sol1:=dsolve({eq,y(0)=0}, numeric, parameters=P):
sol1(parameters=[1,2,3,4]):
dsolve([eq, y(0) = 0], y(x), type = 'series'):
sol2:=convert(%, polynom):
sol2:=eval(rhs(sol2),P=~[1,2,3,4]);

plots[display](plots[odeplot](sol1,[x,y(x)], x=0..5, color=red, legend=`Numerical solution`), plot(sol2, x=0..5, color=blue, legend=`Solution by series`), view=[0..5,-1..1.2]);

In  #2  should be

convert~(L, string);

                               ["Norman.Mailer", "Richard.Brautigan"]

Even simplier example:

g:= a -> int(x+a, x=a..2*a):

g(1);   # OK

eval(g(x), x=1);  # The error ( a premature calculation)

eval(g(z), z=1);  # OK

eval('g(x)', x=1);  # OK (a workaround)

V:=Int(x*sqrt(2*x^4+3), x);

value(V);  # The result (a direct calculation)

IntegrationTools[Change](V, u = sqrt(2)*x^2);  # Calculation by a change

value(%);

combine(expand(eval(%, u = sqrt(2)*x^2)));   # The result

Int  is an inert form of an integral.

Here is a simple example of the equation with the explicit solution, in which the same error:

dsolve({diff(y(x),x)*(x-1)^2+1=0, y(0)=0});

sol:=dsolve({diff(y(x),x)*(x-1)^2+1=0, y(0)=0}, numeric);

sol(2);

 It is obvious that the curve going from point (0,0) can not be continuously extended further right of  x = 1

plottools:-getdata  command seems appeared in Maple 15. For older versions instead of  plottools:-getdata(A)[3];  you can write

op([1, 1], A);

First, we find 3 specific points on the surface, then - the equation of the plane through these points, and then we prove that each point of the surface lies on this plane:

restart;

r:=[(u-v)^2, u^2-3*v^2, (1/2)*v*(u-2*v)]:

P1, P2, P3:=seq(eval(r,p), p=[[u=0,v=0],[u=1,v=0],[u=0,v=1]]);

LinearAlgebra[Determinant](Matrix([[x,y,z], P2, P3]))=0;

expand(eval(%,[x,y,z]=~r));

Your system is inconsistent for any  rho<>0 . First we eliminate  rho  from equations 2 and 3 and then solve the system of two equations:

restart;

sys:=[8*g(t)^3*diff(g(t),t$2)+4*(g(t)*diff(g(t),t))^2+1=0,rho=-1/g(t)-2*(diff(g(t),t)+t*diff(g(t),t$2))-t/(2*g(t)^3),rho=(-t/g(t))*(diff(g(t),t))^2+t/(4*g(t)^3)];

R:=eliminate(sys[2..3], rho);

Sol:=dsolve([sys[1],op(R[2])]);

eval(sys,Sol[1,1]);

eval(sys,Sol[1,2]);

 For  rho=0  solutions are  Sol .

In the procedure  linalg  package was used but it is not called. Add the line  uses linalg; after local-line of your procedure:

Cocycle:=proc(L,n)

local i,j,k,h,v,u,w,C,eqns,e,f,g;  

uses linalg;

... 

nm

At first I misunderstood your demands. The new procedure  P1  solves the problem for all selected variables in any powers (including the symbolic powers). Note that  simplify  with sideral relation does not work for symbolic powers.

restart;

P1:=proc(expr, S::list)

local Rule;

Rule:=product(S[i]^p[i]::anything, i=1..nops(S))=0;

applyrule(Rule, normal(expand(expr)));

end proc: 

Examples of use:

f0 := a^4+4*a^3*b+6*a^2*b^2+4*a*b^3+b^4:

f1 := 3*(a*b-2*c):

f2 := (a*b-2*c)/(-a*b+c):

f3:=a^b*c:

f4:=a+b*c:

f5:=c+a^m*b^n:

P1(f0,[a,b]); P1(f1,[a,b]); P1(f2,[a,b]); P1(f3,[a,b]); P1(f4,[a,b]); P1(f5,[a,b]);

restart;

P := expr->evalindets(normal(expand(expr)), And(`*`, satisfies(t->nops(indets(t))>1)), t->0): 

Examples of use:

f0 := a^4+4*a^3*b+6*a^2*b^2+4*a*b^3+b^4:

f1 := 3*(a*b-2*c):

f2 := (a*b-2*c)/(-a*b+c):

P(f0);  P(f1);  P(f2);

This works for any rational expressions in several variables, if the denominator is not equal to 0 after substitutions.

In your code you are using a point instead of a multiplication sign many times. For example compare the outputs:

simplify(sin(t).sin(t)+cos(t).cos(t));

simplify(sin(t)*sin(t)+cos(t)*cos(t));

In Maple a point is used for the multiplication of matrices and vectors, e.g.

V:=<a, b, c>:

V.V  assuming real;