@Jjjones98 The first step is to show that if f is any polynomial of degree 2*n-1 or less (n is 4 in your example), then the approximation is in fact exact equality. For further details, see the Wikipedia article.

@Jjjones98 Yes, the approximation relationship that you state holds for all suitably smooth functions f (but not for just any f). However, I don't know how you'd use Maple to show that in general. If you want to use Maple to show it for specific examples of f, that's easy. For a discussion of the general case, I suggest that you read the Wikipedia article "Gaussian quadrature".

Here's an example. I also compute the coefficients c_k, the formula for which I got from the Wikipedia article.

Download Guassian_Quadrature.mw

@ebrahimina You changed x[j] from my code to x[i]! The x[j] comes directly from the formula in your Question, so how could you make this mistake?

Kitonum's code is essentially the same as mine. How is it that you transcribed his correctly, yet you've transcribed mine incorrectly twice? Are you trying to make me look like a fool?

@assma You have

local x, y; z,

which should be

local x, y, z;

@torabi Your equation has diff(u(r),r) raised to to the powers n-1 and n-2 with n < 1. Your boundary conditions say that this derivative can be 0. So, it's dividing by 0.

@Jjjones98 What do you mean by the "given" integral? You didn't "give" any integral in your question.

@gaurav_rs If you're talking about solving multiple systems--- A1.X1 = B1, A2.X2 = B2, etc. ---that can be done over the 24 processors by distributed processing, which is a form of multi-processing that is more loosely structured than parallel processing. Specifically, the concurrent processes don't share memory. It is done by Maple's Grid package. It is much easier to implement than parallel processing because you don't need to control multiple processes trying to write to the same areas of memory. But likewise, it is slower than parallel and uses much more memory.

@gaurav_rs Two hours seems like way more than what should be required, unless there's something important that you're not telling me, such as that you want thousands of digits of accuracy. Here are my results:

Digits:= 30:
n:= 400:
A:= LinearAlgebra:-RandomMatrix(n$2, datatype= float):
B:= LinearAlgebra:-RandomVector(n, datatype= float):
X:= CodeTools:-Usage(LinearAlgebra:-LinearSolve(A, B, method= hybrid)):
memory used=0.66GiB, alloc change=2.45MiB, cpu time=4.28s, real time=4.89s, gc time=1.92s
LinearAlgebra:-Norm(A.X - B, infinity)*n/LinearAlgebra:-Norm(B,1);

1.69177744953225012309207287051*10^(-27)

So, about 27 digits of accuracy in under 5 seconds.

As far as I know, there is no pre-packaged parallel method for LinearSolve; you'd need to write your own.

Many thanks to all who commented on my code, and especially to VV for his extensive testing for accuracy. Here is some improved code. The numeric computation runs in pure evalhf mode, so it's much faster, and it's probably more accurate also. I got 6 correct digits starting at the millionth position, the accuracy verified from the table at the end of the paper. I also simplified the conversion to a hex string to a single sprintf.

It should be pretty straightforward to convert this to Matlab.

SSj:= proc(d::nonnegint, j::posint)
local k;
   Frac(
      Frac(add(ModExp(16, d-k, 8*k+j)/(8*k+j), k= 0..d)) 
      + add(1/(8*k+j)/16^(k-d), k= d+1..d+14)
   )
end proc:

ModExp:= proc(b, N, m)
#Compute b^N mod m
local t:= 1, n:= N, r:= 1;
   while t < n do t:= 2*t end do;
   do
      if n >= t then r:= irem(b*r, m); n:= n-t end if;
      if t=1 then return r end if;
      t:= iquo(t, 2);
      r:= irem(r^2, m)
   end do
end proc:
      
Frac:= proc(x::numeric)  
local f:= x - trunc(x);
   `if`(f >= 0, f, 1+f)
end proc:

Pi_BaileyBorweinPlouffe:= proc(d::posint)
#Prints hexadecimal digits of Pi starting at an arbitrary position, d.
#9 digits will be printed because
#that's the limit of hardware-float arithmetic according to the paper.
#The accuracy of all 9 digits is not guaranteed.
local x, dm1:= d-1;
   Digits:= trunc(evalhf(Digits));
   x:= evalhf(16^9*Frac(4*SSj(dm1,1) - 2*SSj(dm1,4) - SSj(dm1,5) - SSj(dm1,6)));
   sprintf("%09X", round(x))
end proc:

What's the approximate size, n? What's the type of the coefficients (floats, integers, etc.)?

As far as I know, the only thing available for scheduling is Threads:-Sleep.

I'm not familiar with DifferentialGeometry, so I can't answer your question. I just wanted to point out that your posted file shows Maple giving a response to each of your commands. Perhaps those responses are wrong, but that's different than Maple not responding to your commands.

@briceM Sorry about that. The M -~ -1 needs to be changed to M -~ (-1). I already corrected the code in the Answer, so you can just copy-and-paste that.

Considering that you have numerous terms with the fractional exponents 2.41, 1.41, or 0.41, it seems likely that solutions are complex. If you remove the range restrictions of the variables and include the complex option, you'll get solutions.

@Lali_miani The set E is ordered, but the order is chosen by Maple for its own convenience. I'm guessing that your version of Maple is before Maple 12, and thus the order is based on addresses and may vary between sessions. Later versions of Maple use a set ordering that doesn't vary between sessions, but may vary between versions. You should never rely on the set order except to know that sets will be ordered, but in a way that you can't control.

@Kitonum Imagine looking up a word in a dictionary or a name in a phone book. I mean actual books, not something on a computer. Surely you'd make use of the fact that the words or names are sorted rather than looking at every entry. Checking that the entries are sorted would also be a huge waste of time. Rather, we must assume that they are sorted. If our search happens to prove that assumption false, maybe we report that or maybe we fail silently; that's not very important.

Being someone who answers many Questions here, I think that sometimes we should consider why the Question is being asked, what "bigger picture" is the asker trying to learn, what some unseen third party (such as an instructor) hopes the asker will learn, etc. That's what I mean by pedagogic spirit. After all, don't we teach the same courses as those unseen instructors?