I can't see exactly what is causing that error message without seeing your Phi---I can't duplicate the message---but try 

int~(Phi^%T*f, x= 0..1);

Note that ScalarMultiply is unnecessary: The operation is handled by ordinary `*` multiplication.

The tilde (~) after the int is needed to map the operation over a vector.

You could use mu[0,1] or mu[`01`] or mu[o1] or mu[O1].

Use readline to read the first line. Use ImportMatrix to read the rest. ImportMatrix has an option to skip the first n lines of the file.

My experience is that Grid operations do not work very well when the operation being replicated is trivial. In this case, the symbolic application of ln, the operation is utterly trivial, just an unevalated fuction call. My guess is that each step should take at least 1/10 second to use Grid.

Try adding the option linestyle= solid to the plot command.

f''(0) is what your last line of code prints. It's what you called X2.

Most random values are 12 digits. Try

eval(eqn1, y= rand());

You will almost certainly get a numeric overflow. I believe that testeq is somewhat naively doing this. The safer way would be to use &^ instead of ^. I believe that testeq is using modular arithmetic with a large random prime modulus, but it is not changing ^ to the safer &^.

The above is just my speculation based on reading ?testeq.

You can avoid the cut-and-paste step: In file 1.mw, after the definition of aver, include the command

save aver "aver.mpl";

You can use any filename that you want in the quotes. Make sure that the command is executed; this won't work if the command is just sitting there unexecuted. In file 2.mw, use the command 

read "aver.mpl":

This part is the same as TomLeslie's Answer. (The parentheses in his read command are superfluous.) Indeed, the save command creates the plain text file for you.

int(ln(1+sqrt((m+1)/m)), m);
simplify(%, symbolic);
value(%);

While the question may seem silly as posed, it is not silly to count (either exactly or by estimate) the solutions. Generalizing, it is very interesting to count the integer lattice points in a bounded convex high-dimension polytope. For the above problem, there are about 14 billion solutions. The exact number is 13,971,037,148. I computed this with the code below in 20 minutes real time on an i7 with eight CPUs, using parallel processing.

It was quite tricky to make this code efficient. Amazingly, many of the fundamental functions for integer programming are not threadsafe: abs, ceil, floor, round, frac. The function frem is threadsafe, but not evalhfable. That leaves only irem, iquo, and trunc as safe to use. So, I had to write a replacement for ceil using these. My procedure Ceil below is called inside the most deeply nested loop. A very small tweak in its efficiency would have a large effect on the overall efficiency. Also, it would help if it could somehow be compiled.

It is amazing that the function ceil is still not threadsafe when used inside evalhf code! I verified this by trying it and getting the wrong count of solutions. How can that be? Isn't the evalhf code supposed to be using a separate copy of ceil that avoids all symbolics?

Another way to improve the efficiency would be to unroll the recursion so that evalhf could be applied at a higher level. Doing this by hand makes the code too ugly to be bearable. But perhaps the unrolling could be done automatically via some metaprogramming.

All comments welcome.

Download LatticeCount.mw

You need to take a few steps to complete this integration:

1. Don't use the same variable as both the integration variable and as a limit of integration. I changed your upper limit of integration from x to X.

2. Assume that X is in some interval. I chose [0,5] somewhat arbitrarily.

3. Use expand to factor the functions of t outside of the integral.

In the code below, I used D notation for your integrand simply because it's less to type. You'll get the same answer if you use the diff notation that you had.

assume(X > 0, X < 5);
int(D[1,2](w)(x,t)^2 + D[1](w)(x,t)*D[1,2](w)(x,t), x= 0..X);
expand(%);

It is not a bug. Your pp has a root at cp=3220 for any value of x. You can see this by eval(pp, cp= 3220), which is identically 0. Because of floating point variation, Roots will find that root at slightly less than 3220 or exactly 3220 or slightly more than 3220. Often Roots will not find this root at all, because the curve becomes vertical (is not differentiable) at 3220. Because you asked for roots in the interval 1..3220, the root will be reported when it is less than or equal than 3220. Since you are not interested in seeing the root at 3220 on your graph, change the interval to 1..3219 and you will never see them.

Edit: Added sentence about Roots often not finding the root at 3220.

1. You seem to have some confusion as to what "leading coefficient" means. It means the coefficient of the term with the largest exponent on the main variable. What you are calling the leading coefficient is actually the "trailing coefficient."

2. Polynomials can be sorted into ascending or descending powers of their main variable with the sort command. For example:

sort(collect(expand(lhs(e4)),t), t, ascending);

3. You cannot control the order that items appear in a Maple set. If you use curly braces {}, it will be a set. If want to maintain an order, use square brackets [ ] to make it a list.

4. Why do you care about the order? Relying on a certain order is often (but not always) a sign of poor coding practice.

Use convert to convert to diff form and another convert to convert back to D form. It can all be done in one line.

eq:= D[2](u)(x,t)= a^2*D[1,1](u)(x,t):
expr:= D[2,2](u)(x,t):
convert(eval[recurse](convert([expr, [eq]], diff)[]), D);

Your immediate problem is that you can't use a formal parameter, x, as the index of a for loop. I think that you meant for i to be the index.

There are other problems which don't produce syntax errors: By local top::0, you surely mean to initialize top to 0. That should be top:= 0. And by 

top =top+(((x^(2+i))*top)^(1/(n+2-i)))

you surely mean to update the value of top. So that should be 

top:= top+(x^(2+i)*top)^(1/(n+2-i))

Now the problem is that the above line can never change top from 0. So your formula is wrong, but I have no idea how to correct it because I don't recognize what this formula is supposed to do.

With all my corrections, the code becomes:

soru:= proc(n, x)
local i, top:= 0;
     for i from 2 to n do
          top:= top+(x^(2+i)*top)^(1/(n+2-i));
          print(top);
     od;
     top
end proc;