Maple distinguishes between a matrix with m rows and 1 column and a column vector with m elements.

You can convert the matrix into a vector, which VectorCalculus:-Jacobian will accept.

Try this:

v:=<x^2*y,y+x>;
M:=convert(v,Matrix);
type(M,Matrix);
type(v,Matrix);
VectorCalculus:-Jacobian(v,[x,y]);
VectorCalculus:-Jacobian(M,[x,y]);
w:=convert(M,Vector);
VectorCalculus:-Jacobian(w,[x,y]);

You can show the two (or more) animations together with display.

An example (not very exciting):

with(plots):
p1:=animate(plot,[sin(t),t=0..T],T=0..6*Pi):
p2:=animate(plot,[cos(t),t=0..T,color=blue],T=0..6*Pi):
display(p1,p2);

If you use assumptions on y, you are OK:

restart;
g := (1-x) * ln(sqrt(x^2 + y^2 + 1) +1 );
h:= subs(y = 0.5, g);

h1 := int(h, x = 0..1);

g2 := int(g, x = 0..1) assuming y>0;
h2 := evalf(subs(y = 0.5, g2));
simplify(h2);

assuming y::real works too, but the result is more messy.

Maple is not doing numerical integration in the session above, not even if wrapped in an evalf. Try setting infolevel[int]:=1: before doing the first integration. If you want numerical integration in the first case, you should do

h1 := evalf(Int(h, x = 0..1));

Notice 'Int' instead of 'int'.

M:=Matrix([[1,1,0],[0,0,1],[1,0,1]]);
Ind:=Matrix(3,(i,j)->if M[i,j]=1 then [i,j] else NULL end if);
convert(Ind,list);

After correcting the syntax in your initial conditions you still have a problem.

You have two equations determining diff(psi(t),t,t), eq1 and eq3. Together with eq2 they seem to be inconsistent.

restart;
eq1:=(l^2+r^2/2)*(sin(theta(t))^2*diff(diff(psi(t),t),t)+2*diff(psi(t),t)*diff(theta(t),t)*cos(theta(t)))-r^2*(cos(theta(t))*diff(psi(t),t)+omega(t))*diff(theta(t),t)=0;
eq2:=(l^2+r^2/2)*(diff(diff(theta(t),t),t)-(diff(psi(t),t))^2*sin(theta(t))*cos(theta(t)))+r^2*(diff(psi(t),t)*cos(theta(t))+omega(t))*diff(psi(t),t)*sin(theta(t))=g*l*sin(theta(t));
eq3:=diff(diff(psi(t),t)*cos(theta(t))+omega,t)=0;
l:=0.05;r:=0.063/2;g:=0.81;omega:=evalf(2*Pi*10000/3600);
ic1:=theta(0)=evalf(9*Pi/180),psi(0)=0,D(theta)(0)=0,D(psi)(0)=0;

eq1;eq2;eq3;

#This doesn't succeed:
Sol:=dsolve({eq1,eq2,eq3} union {ic1},{theta(t),psi(t)},numeric);
Maybe there is a good reason. From eq1,eq2, and eq3 you can find two odes, one of which is second order in theta (L2) and the other (L13) is first order in psi. These two can be solved: 
L1:=op(expand(solve(eq1,{diff(psi(t), t, t)})));
L3:=op(solve(eq3,{diff(psi(t), t, t)}));
L2:=op(solve(eq2,{diff(theta(t), t, t)}));
eval(L1,L3);
(rhs-lhs)(%)=0;
L13:=op(solve(%,{diff(psi(t), t)}));
ic1;
ic2:=ic1[1..3];
Sol:=dsolve({L13,L2} union {ic2},numeric,output=listprocedure);
#However, it appears that the second derivative of psi found from Sol by fdiff and the same found from L3 using Sol don't agree:
TH,dTH,PSI:=op(subs(Sol,[theta(t),diff(theta(t),t),psi(t)]));
fdiff(PSI,[1,1],[1.2345]);
psi2:=eval(rhs(L3),L13);
eval(psi2,{theta(t)=TH(t),diff(theta(t),t)=dTH(t)});
eval(%,t=1.2345);

Your f is not defined as a function, so should not be used as such.

Replace f(x) by f:

s := solve(y = f, x);

#Maple reports 9 solutions:
nops([s]);

#one of which contains no obviously imaginary parts:

remove(has,[s],I);

Since DEtools is not a module you cannot use 'uses' and you should not use 'with' within procedures.

type(DEtools,`module`);
type(DEtools,table);

DEtools is a table (used for packages before modules were introduced in Maple).
Use the long form DEtools[DEplot] inside procedures.

Example.

eq:=diff(y(t),t)=y(t);

p:=proc(t0::realcons,y0::realcons,T::realcons) global eq,y,t;    
    DEtools[DEplot](eq,y(t),t=t0..T,[y(t0)=y0],_rest)
end proc;

p(0,1,2,linecolor=blue);

By experience I know that one has to be extremely careful when using an indexed name like 'a[12]' and also its unindexed form 'a'. However, when they are not assigned values, it should not be as dangerous.

Nevertheless, to avoid speculating about that I replaced the bare form 'a' with 'b'.

I ran the following, and found solutions in terms of RootOf's, which is not surprising. But I only found valid solutions when first converting tanh to exp.

restart;
sys:=[-xt[1,1]+b*tanh(xt[1,1])+a[12]*tanh(xt[2,1])=0,-xt[2,1]+a[21]*tanh(xt[1,1])+b*tanh(xt[2,1])=0];
solve(sys,[xt[1,1],xt[2,1]]);
allvalues(%);
#Solutions false:
eval(sys,{xt[1,1]=ln(I),xt[2,1]=ln(I)});
convert(sys,exp);
res:=solve(%,[xt[1,1],xt[2,1]]);
indets(res,name);
#Testing the solution for concrete randomly chosen parameters in the interval (0,10)
#and looking for nontrivial solutions only.
rnd:=10*RandomTools:-MersenneTwister:-GenerateFloat:
RES0:=[[xt[1,1]=0.,xt[2,1]=0.]]: RES:=RES0:
for k from 1 to 100 while RES=RES0 do
  val:=[b=rnd(),a[12]=rnd(),a[21]=rnd()];
  RES:=simplify(fnormal(evalf(eval(res,{val[1],val[2],val[3]}))))
end do:
RES; val;
k;
eval(sys,{op([1,1..2],RES),op(val)});

To write g_{i,j} is certainly not allowed in 1d-input, which I always use.

But if you do

p:=(k,j)->if k=0 then max(k-j*S,0) else g(k,j) end if;

after having defined the function g in the proper manner (*), then p(1,4) should execute as expected.

(I changed K to k assuming that was what you meant).

(*) Example. g:= (i,j) -> i^2-j*i;

The derivative exists and is 1/2.

A:=Int(cos(1/sqrt(u))^2,u=0..x);

IntegrationTools:-Change(A,u=t^2);

B:=combine(convert(%,radical));

#We need only consider

C:=Int(t*cos(2/t), t = 0 .. sqrt(x));

IntegrationTools:-Change(C,t=sqrt(x)/u,[u]) assuming x>0;

#Now we shall show that
E:=Int(cos(n*u)/u^3, u = 1 .. infinity);

#tends to zero as n tends to infinity.

#We split the integral in two, where M is to be chosen later
IntegrationTools:-Split(E,M);

# and use an integration by parts on the first term:
applyop(IntegrationTools:-Parts,1,%,u^(-3));

The remaining details follow by a standard argument.

There are a few errors in your statement of the problem. sin*t should be sin(t), 'pi' should be 'Pi', u(0,t) should not be assigned to. 

These errors are corrected below.

PDE:=diff(u(x,t),t)+2*diff(u(x,t),x)=0;
bc:=u(0,t)=piecewise(0<t and t<=2*Pi, sin(t),0); 
pdsolve({PDE,bc});
U:=rhs(%);
plot(eval(U,t=0),x=-5*Pi..Pi,thickness=2,caption="The solution for t = 0");
plots:-animate(plot,[U,x=-5*Pi..7*Pi,thickness=2],t=0..3*Pi);

Using the polar form for t, t = r*exp(I*v), we immediately see that v = Pi/5 + p*2*Pi/5, where p = 0, 1, 2, 3, 4.

Then

restart;
eq:=(-t)^5=abs((t^2+t+1)/(-t));
evalc(eval(eq,t=r*exp(I*(Pi/5+p*2*Pi/5)))):
eq3:=combine(%) assuming p::integer,r>0;
#Trying p=0:
solve({eval(eq3,p=0),r>0},r);
evalf(%);
#Now trying all p's:
[seq(solve({eq3,r>0},r),p=0..4)]:
L:=evalf(%);
nops(L);
L2:=subs~(L,r);
#Maple found nothing for p = 2, which corresponds to t = -1.
#We consider p = 2 separately:
solve({eval(eq3,p=2),r>0},r);
factor(eval(eq3,p=2));
simplify(%) assuming r>0;
solve({%,r>0},r);
res:=zip(`*`,L2,[seq(exp(I*(Pi/5+p*2*Pi/5)),p=[0,1,3,4])]);
plots:-complexplot([-1,op(res)],style=point,symbolsize=20);

It seems from Maple that the only solution is -1.

I did the following:

eq:=(-t)^5=abs((t^2+t+1)/(-t));
res:=solve(eq,t);
op(0,res);
evalf(res);
#None of the following are satisfied for obvious reasons!
map(evalf,op(1,res));

#The seven values:
evalf(op(2,res))

The answer from solve is certainly not perfect. It comes in the form of a piecewise expression.

Maple finds 7 t-values (6 RootOf's and -1). Maple stalls on the condition that the fifth powers of the RootOf's be negative, which they aren't: They are not even real.

Apart from the stalling, if Maple is right, then the only solution is -1.

The problem is that y is used in two different ways. y is first implicitly defined as a table by the assignment y[1]:=0.5:

Later it is used as just a name in

s:=solve(y=g[i](x),x):
h[i]:=unapply(s[1],y);

If you change 'y' to e.g. 'yy' the loop runs.

The use of a table name as a mere name is likely to cause strange problems. The resulting error messages are usually very opaque.

If you could provide us with the actual system and initial conditions then we wouldn't have to make guesses.