LinearAlgebra:-Eigenvalues is not in the list at ?threadsafe, so there's no point trying to use it in multithreaded code.

Change rho[0]=1.33 to rho[0]:= 1.33. Likewise for T[0].

A search for "quaternions" at the Maple Applications Center turns up a handful of packages.

This type of problem can be done with my LogicProblem package, available at the Maple Applications Center, although in some ways it feels like using an elephant gun to swat a fly. One advantage over the other methods is that if the solution is unique, then my code will prove it unique.

Download TopoSort.mw

In order to connect the points, they must be sorted. The code below does that. The centroid of each group is found, then angle that each line segment from point to centroid makes with the positive x-axis is used to sort.

Download StarPlot.mw

First, you have to decide whether you want to measure the real (aka elapsed, clock) time or the processor (aka execution, cpu) time. In the above code, I think it makes the most sense to measure the processor time, which means that you should use time() rather than time[real]().

Second, you have to realize that the resolution of the time() command is 1/64 of a second, which is reported as 0.016. So, for a timing to be meaningful, it needs to be at least several times this value.

Third, you need to account for the effect of garbage collection on your timing. If a garbage collection is only done once during your timing, that garbage collection may be taking a significant percentage of the time. In order for the timing to be meaningful, garbage collection needs to happens either 0 times or several times during the timing. You can count the garbage collections with kernelopts(gctimes):

gc(): # Force collection
stgc:= kernelopts(gctimes):
st:= time():

     ... test code here

gct:= kernelopts(gctimes):
if gct <> stgc then  gc()  end if:
time() - st;

Fourth, there may be a "first-time effect". The first time the code is run (first time after a restart), there may be time involved in reading the Maple procedures from libraries. This time may not apply during subsequent runs performed without doing an additional restart. Also, memory needs to be allocated for the first time. Therefore, you need to compare first-time timings or subsequent-time timings; you can't compare a first-time timing with a subsequent-time timing.

Fifth, if you really do want to measure real time, then your code---including the timing commands---must all be in one execution group. Otherwise, you're measuring the time it takes you to press the keys, which is irrelevant.

Finally, despite my best efforts to eliminate confounding factors, I still often get different timings for the same piece of code. There are confounding factors that I don't understand or don't know how to control. You may need to take a modal value from several timings. Paying attention to the above factors will at least reduce the variation.

If you want to model a spring/damper system with variable damping as a differential equation, then you can change the damping based on the velocity. I don't know how you can change it based on acceleration. Here's an example:

Download damping.mw

The operation can be done entriely within one Matrix, that Matrix representing your spreadsheet. Let's say your data is in the 3rd to 5th rows, and you want to multiply columnns 1 and 2 and store the result in column 3:

A:= LinearAlgebra:-RandomMatrix(5, generator= rand(-9..9));

A[3..5, 3]:= A[3..5, 1] *~ A[3..5, 2]:
A;

As it stands, you cannot guarantee that the values of y used in the 28x2 are a subset of the values of y used in the larger two matrices. So trimming the larger two matrices down to the size of the smaller may not help. The way around this is to use option adaptive= false to the plot command. Then the three matrices will be the same size and have the same second column (the y column).

Here's an example of making it into a single four-column matrix:

plot(
     [[y,y,y=0..3], [y^2,y,y=0..3], [y^3,y,y=0..3]],
     numpoints= 25, adaptive= false
);

data:= [plottools:-getdata(%)];

< op([1,-1], data)[.., [2,1]] |
    op([2,-1], data)[.., 1]       |
    op([3,-1], data)[.., 1]
>;

Round:= (x,n)-> parse(sprintf(sprintf("%%.%df",n),x));

where x is the value to be rounded, and n is the number of decimal places. For example Round(350.4971, 2);

It's as easy as

plot(A[.., [1,3]]);

The .. means "all rows".

The help page ?dsolve,numeric,bvp just says that the midpoint methods can help with "harmless endpoint singularities". Here's an an example: y'' = sin(x)/x on 0..1. The derivative expression cannot be properly evaluated at the left endpoint because of the division by zero, but we know that this singularity is removable because limit(sin(x)/x, x= 0) = 1. Since the midpoint methods evaluate the derivative expression at the midpoints of the subintervals rather than the endpoints, they will never use x=0 in this problem---0 will never be a midpoint of a subinterval.

Download midrich.mw

That's right: You cannot integrate numerically when the integral has a parameter, the kk in your case. The integral will be done if you specify a value for kk. Considering the complexity and non-smoothness of the integrand, you need to increase your error tolerance from the default in order to get an answer in reasonable time. This is done with the epsilon option.

JJJ:= Int(max(0., (0.9483573506e-3*(-1.*sin(a)*cos(w)-1.*cos(a)*sin(w)*sin(b)))*cos(a)*cos(b)^2*(-58.5*signum(cos(b)*sin(w)*sin(b))*kk+200.*cos(b)*sin(w)*sin(b))*Heaviside(-58.5+200.*cos(b)*sin(w)*sin(b)*signum(cos(b)*sin(w)*sin(b))*kk)*Heaviside(1.-.98*cos(b)^2)/sqrt(1.-.98*cos(b)^2)) +
min(0., (0.9483573506e-3*(-1.*sin(a)*cos(w)-1.*cos(a)*sin(w)*sin(b)))*cos(a)*cos(b)^2*(-58.5*signum(cos(b)*sin(w)*sin(b))*kk+200.*cos(b)*sin(w)*sin(b))*Heaviside(-58.5+200.*cos(b)*sin(w)*sin(b)*signum(cos(b)*sin(w)*sin(b))*kk)*Heaviside(1.-.98*cos(b)^2)/sqrt(1.-.98*cos(b)^2)), [a= 0..2*Pi, b= 0..Pi/2, w= 0..2*Pi], epsilon= 1e-4, method= _cuhre):

evalf(eval(JJJ, kk=1));

This answer should only be trusted to 4 significant digits.

Using unprotect let's you assign to D; it doesn't change or remove the value that D already has. To do that, you need unassign in addition to unprotect:

save_D:= eval(D):  #Just in case you need it later.
unprotect(D):
unassign('D'):

However, if you have Maple 17, it would be better to use the local option (although not as described by Georgios):

local D;  Lett:= D;

This leaves the original D accessible as :-D.

Here is your program translated into modern Maple.

Download Matrix_Calculation.mw