When you don't have a suitable t variable to use with animate or if your frames are already created, then you use display with the insequence option to produce the animation. Simply change your last line to

plots:-display([seq(frame[j] $ iquo(50, lastj+1), j= 0..lastj)], insequence);

To save the animation as a GIF, right click on the animation, and select Export -> GIF, which is what I did to copy it below.

I think that your limits of integration are best handled in spherical coordinates. I will exchange the roles of y and z in your integral because it is easier for me to think of the hemisphere z >= 0 than y >= 0. The integrand z becomes rho*cos(phi) and the differential dx*dy*dz becomes rho^2*sin(phi)*dphi*dtheta*drho.

int(int(int(rho*cos(phi)*rho^2*sin(phi), phi= 0..Pi/2), theta= 0..2*Pi), rho= 0..1);

You're opening a can of worms by claiming that there is a "first" die. If you truly mean that, then the answer is

S:= {seq([5,j], j= 1..6)}:
E:= select(r-> `+`(r[]) > 9, S):
nops(E)/nops(S);

If what you actually mean is "given that at least one of the dice is a 5", then the answer is 

S:= {seq([5,j], j= 1..6), seq([i,5], i= 1..6)}:
E:= select(r-> `+`(r[]) > 9, S):
nops(E)/nops(S);

restart:

Curve:= t-> [a*cos(t), a*sin(t), b*t]:
Surface[xyz]:= x^2 + y^2 = a^2:
Surface[parametric]:= (s,t)-> [a*cos(t), a*sin(t), s]:

#Is the curve on the surface?
is(eval(Surface[xyz], [x,y,z]=~ Curve(t)));
                                          true

#Is the parameterized surface the same as the surface in xyz form?
is(eval(Surface[xyz], [x,y,z]=~ Surface[parametric](s,t)));

                                true

a:= 1:  b:= 1:
plots:-display(
     [plot3d(Surface[parametric](s,t), t= 0..2*Pi, s= 0..2*Pi*b),
      plots:-spacecurve(Curve(t), t= 0..2*Pi, thickness= 4, color= red)
     ],
     axes= none, orientation= [10, 31, 22]
);

If you end a command with a colon, it suppresses the output of the command. To see the output, use a semicolon.

If, for example, you want to do a one-tailed test with 12 degrees of freedom and a 95% confidence interval, then you get the critical value of t by using

Statistics:-Quantile(StudentT(12), .95);

for a two-tailed test, change the second parameter to .975 (since 1 - (1-.95)/2 = .975).

To get the p-value of the test, compute your t statistic, t, then

p:= 1 - Statistics:-CDF(StudentT(12), t);

In order to provide a more detailed answer, such as how to compute the t statistic, I need more details of your problem. What do you mean by "dependent means"? Do you mean that the samples are paired?

I usually prefer eval over the methods described by Thomas Richard:

dd := fsolve({eq1, eq2}, {x[1], x[2]});
evalf(eval(sin(x[1]), dd));

This lets you keep x[1] and x[2] as symbolic variables.

This is a job for fsolve, not solve. You should've uploaded the entire system so that I could test it. You may also try DirectSearch, which needs to be downloaded from the Maple Applications Center.

If you had only allowed non-strict inequalities, then this would simply be a matter of calling Optimization:-Minimize with a trivial objective function. Allowing strict inequalities, it's a little more complicated. Here's a procedure for it:

FeasiblePoint:= proc(
     Cons:= set({`<`, `<=`, `=`}),
     {epsilon::positive:= 10.^(2-Digits), epsilonmax::positive:= 1}
)
local
     newCons:= subsindets(Cons, `<`, ineq-> lhs(ineq) <= rhs(ineq) - epsilon),
     Sol:= Optimization:-Minimize(0, newCons)[2]
;
     if not hastype(Cons, float) then Sol:= convert(Sol, rational) end if;
     if andmap(is, eval(Cons, Sol)) then Sol
     elif epsilon < epsilonmax then thisproc(Cons, 'epsilon'= 10*epsilon)
     else error "No feasible point found; epsilonmax exceeded."
     end if
end proc:   

Example of use:

FeasiblePoint({ x < 1 - y, y < 0, y + 1 < x });

If kx=0 or ky=0 or both, then the integral is 0, as can be determined by direct evaluation. The result returned for the general case is still valid for kx=0 and/or ky=0. So there's no need for piecewise.

The absolute value below is in case you don't know which function is on top. That's trivial in this case, but I wanted to be a bit more general.

F:= (-x^2+9) - (x+3):
abs(int(F, x= `..`(solve(F, x))));

TangentLine:= proc(f::algebraic, at::name= algebraic)
local x, a;
     (x,a):= op(at);
     eval(diff(f,x), x= a)*(x-a) + eval(f,x= a)
end proc:

NormalLine:= proc(f::algebraic, at::name= algebraic)
local x, a;
     (x,a):= op(at);
     eval(f, x= a) - (x-a)/eval(diff(f,x), x= a)
end proc:

f:= ln(3*x)+3; a:= exp(1);
                        f := ln(3 x) + 3
                          a := exp(1)

TL:= TangentLine(f, x= a);

NL:= NormalLine(f, x= a);

plot([f, TL, NL], x= 1/2..10, y= 1/2..10, legend= [f, tangent, normal]);

The following computation strongly suggests that the value that you seek is undefined:

 MG:= MeijerG([[-.3+eps], []], [[.8, 1.3, -.8, -.3, -1.3], []], 1.):
eval(MG, eps= 10^k) $ k= -6..-2;

eval(MG, eps= -10^k) $ k= -6..-2;

Like this:

Ys:= [1,2]: #List all y values here.
Zs:= [3,4]: #List all z values here.

G:= (y,z)-> sqrt(D[1](x)(y,z)^2*dy^2 + D[2](x)(y,z)^2*dz^2):
zip(G, Ys, Zs);

See the command ?map. Example:

V:= <1,2>;

map(tan, V);