@hopereza 

Hmm, the tilde is not working in your worksheet. I don't understand why yet. But I tried this in your worksheet, and it definitely works:

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

@Markiyan Hirnyk I converted your Reply to an Answer and gave it a vote up. I think that your solution is better than mine.

@Rouben Rostamian  ScalarMultiply multiplies a Matrix or Vector (its first argument) by a scalar (its second argument) according to ?ScalarMultiply.

@itsme In Maple 18, I find no command named Grid:-Set. There is a Grid:-Setup, but this does not seem to have anything to do with copying memory.

@smith_alpha 

C:= ImportMatrix(
     "C:/Users/Carl/Desktop/Matrix.txt",
     skiplines= 1, source= delimited, delimiter= " ", transpose
);

@mthkvv Threads has no problem with trivial operations. In the case of Grid, the entire memory space must be copied to each process, plus a new kernel has to be created.

@smith_alpha 

Read the first line with 

Z:= fscanf(filename, "%d")[];

Technically, the whole file can be read with fscanf because there are scan codes for matrices. But it'll be easier to use ImportMatrix, which just constructs the correct fscanf command for you.

@ 

You wrote:

If Maple solves with rkf45dae it probaly means that is solves dy/dt =f with algebraic equations may be as a constraint.

That is not the understanding that I get from the help page ?dsolve,numeric,DAE_extension. In particular, it explicitly says that the systems can be partially implicit.

Patient, convert all of your equations to dy/dt =f .

That doesn't appear to be possible for this system.

If needed differentiate the algebraic equations,

Patient's system doesn't have any algebraic equations. Nonetheless, it seems to require a DAE solver because of its implicit nature. It is possible that another of Maple's IVP solvers can handle it. I haven't tried yet.

If some one posts *.mws worksheet that works, I will try. ( or at least in that format, 1D math)

The worksheets posted so far in this thread have been in 1D form.

@acer 

Thanks, Acer. Using option threadsafe and using eval to run compiled code inside evalhf were the things that I didn't know about. I made the outer, recursive procedure evalhfable with your suggestion of passing the arguments in an Array and I declared most arguments of the inner, compiled procedure as float. (This was much much faster than the default integer arguments.) With these changes, I made a factor-of-three improvement of your version, coming in at just under six minutes. And I got 100% processor utilization (M18/64 on Windows 8/64). I am curious if you can get 100% utilization on Linux with this. (I also had a version where I only got 65% utilization.) The worksheet is attached.

Note that I was able to use ceil in the compiled code without any threadsafe issues. It makes sense since the compiled C code should be using a completely different ceil. I'd still like to know how evalhf has threadsafe issues when using ceil.

One thing that was a disaster was trying to have one compiled procedure call another compile procedure. In particular, I compiled a version of Ceil, and ran it in parallel. I used threadsafe in both (and they truly were threadsafe). The execution deadlocked immediately: No processor usage, no memory comsumption, no way to kill the program other than killing the kernel. In another case, I didn't use threadsafe in either. This code executed correctly, but much slower than if I had just used uncompiled Maple. So, do you know anything about one compiled procedure calling another?

Download LatticeCountmodif4f.mw

My algorithm is incorrect. (It has a fundamental flaw, not just an error in coding.) The one posted by Kitonum in a separate Post is correct and much faster. So this is probably a moot point. Still, I learned a lot about making programs faster and was able to combine three of the major techniques (evalhf, Compile, and Threads) into one short program.

@momoklala 

One problem is the B in Eq1. It should be B*. The way it is now makes it a function rather than a coefficient. However, making this change, I still don't get what your instructor got. You may still be entering parts of Eq1 incorrectly. Perhaps you should post the original statement of the problem.

@acer 

Okay, I see that the problem is that eval may not recognize expressions that become 0/0, and it may evaluate them based just on the numerator. So what can be done in this specific case? My unconfirmed belief is that curve crosses the axis at 3220 and its tangent is vertical there. Is there any way to confirm or disprove? 

@ecterrab But PDEtools:-dsubs doesn't work with expressions in D notation, which is the point of this Question.

@Markiyan Hirnyk That 3220 is a root for any x can be verified by exact computation. I explained that in my Answer.

pp:= -55471918776960000*tanh((1/3220)*sqrt(10368400-cp^2)*Pi*x/cp) +
           5350094400*tanh((1/3220)*sqrt(10368400-cp^2)*Pi*x/cp)*cp^2 -
           129*tanh((1/3220)*sqrt(10368400-cp^2)*Pi*x/cp)*cp^4 +
           2670899840*tanh((1/6450)*sqrt(41602500-cp^2)*Pi*x/cp) *
           sqrt(41602500-cp^2)*sqrt(10368400-cp^2):
eval(pp, cp= 3220);

     0

So, it is unquestionable that 3220 is a root. You pollute the computation by using decimal points in the 8000 and 3220. That Roots is unable to find the root at 3220 is up to floating point error. Note that the curve becomes vertical at that point, which can confuse root finders. Check the derivative at 3220.

In Maple 18, I get
Digits:= 150:
eval(x= 8000.0, pp), cp= 3220.0);

     0.

which differs from what you got.

@Sagar With respect to which variable are you integrating from 0 to 1? with respect to t? with respect to X?

@Kitonum Shouldn't you be using chebpade instead of pade?