You almost have it. Your problem is that there is no need for any looping command. Some local variables will help also.

First_Principle:= proc(f::appliable)
local x, h, A:= simplify((f(x+h)-f(x))/h);
     Limit(A, h = 0) = limit(A, h = 0)
end proc:

Example of use:

First_Principle(sin);

You can remove from the solutions those for which the evaluation of R leads to error like this:

Error:= proc() error NumericException(division_by_zero) end proc:
NumericEventHandler(division_by_zero= Error):
CheckSolution:= proc(sol::set(name=anything))
     try
          eval(R, sol)
     catch "numeric exception: division":
          return true
     end try;
     false
end proc:

remove(CheckSolution, sols);

Nothing special is required because of the change of version. Just open it as you would open any other file.

You can discard randomly generated values that don't work, like this

restart:
M:= Array([[a+2*b,b],[x1+y1,z1+z2]]);
pars:= indets(M);
npars:= numelems(pars);
ntrials:= 20:
myRand:= RandomTools:-Generate(list(float(range= -3..3), npars), makeproc):

to ntrials do
     do
          set1:= pars=~ myRand();
          numM:= eval(M, set1);
          if andmap(x-> 0 < x and x < 1, numM) then break end if
     end do;
     LinearAlgebra:-Rank(numM)
end do;

The above takes an average of about 1600 random trials to find each matrix that works.

Apparently the KroneckerProduct command does not work with matrices that have index functions. You need to convert your shape= Circulant matrices to shape= rectangular like this:

A1:= Matrix(A, shape= rectangular):
B1:= Matrix(B, shape= rectangular):
KroneckerProduct(A1,B1);

You may be able to get what you want from ToInert. For example,

ToInert(A = B*C);

odeplot is in the package plots, so you need to refer to it as plots:-odeplot. Or, you load the package with

with(plots):

Then you can refer to the commands in the package without using the plots:- prefix.

PDEtools:-dchange(
     {x(t)= r(t)*cos(theta(t)), y(t)= r(t)*sin(theta(t))},
     {
          diff(x(t),t) = x(t)*y(t)-x(t)^2*y(t)+y(t)^3,
          diff(y(t),t) = y(t)^2+x(t)^3-x(t)*y(t)^2
     },
     {r(t), theta(t)}
);

solve(%, {diff(r(t),t), diff(theta(t),t)});

simplify(%);

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.