Using the package Gravitation, calculating the Einstein tensor (the left hand side of the Einstein field equations) can easily be done as follows:

• Load the package:

  with(Gravitation):
  

• Define/create the metric:

  metricComps := Array(1..4,1..4):
  metricComps[1,1] := +c^2*exp(nu(r)):
  metricComps[2,2] := -exp(lambda(r)):
  metricComps[3,3] := -r^2:
  metricComps[4,4] := -r^2*sin(theta)^2:
  metric := createMetric("LL",metricComps,[t,r,theta,phi]):
  

• Calculate the Einstein tensor, and show its components:

  Einstein := deriveEinstein(deriveChristoffel(metric)):
  Einstein:-getComps();
  

PS: I use my own package because I think that the implementation done in the Maple package tensor is rather confusing because of its non-consistent notation (lowercase letters and uppercase letters with no apparent system), and, especially, because it often requires (as is the case here with the Einstein tensor) one to undertake calculations of many intermediate quantities in order to obtain the final result.

Maybe the following example is useful:

result := Matrix(10,2):
for i from 1 to 10 do
	M := Matrix([[1,i],[-i,1]]);
	result[i,1..2] := LinearAlgebra:-Eigenvalues(M);
end do:
ExportMatrix("... some file path ...",result):

Here is a method that works for any algebraic expression of any number of variables:

getConst := (expr::algebraic,vars::list(name)) -> select(is,expr,freeof(vars)):

where vars is a list of the considered variables of the algebraic expression. An example:

expr := x^2 + y^2 + z^2 - cos(x) - cos(y) - cos(z) + a:
getConst(expr,[x]);
getConst(expr,[y]);
getConst(expr,[z]);
getConst(expr,[x,y]);
getConst(expr,[y,z]);
getConst(expr,[z,x]);
getConst(expr,[x,y,z]);
                  2    2                      
                 y  + z  - cos(y) - cos(z) + a
                  2    2                      
                 x  + z  - cos(x) - cos(z) + a
                  2    2                      
                 x  + y  - cos(x) - cos(y) + a
                         2             
                        z  - cos(z) + a
                         2             
                        x  - cos(x) + a
                         2             
                        y  - cos(y) + a
                               a

It seems that the characteristic polynomial of degree five cannot be solved by radicals. If it could, then (probably) the function allvalues could be used, as given in the following example:

allvalues(RootOf(_Z^3 + 1));

-1, 1/2+1/2*I*3^(1/2), 1/2-1/2*I*3^(1/2) Resorting to numerical methods, maybe the following could be useful:

points := {}:
Fmin := 0:
Fmax := 0.1:
Fnum := 50:
det := Determinant(CharacteristicMatrix(H,lambda)):
for x from Fmin to Fmax by (Fmax - Fmin) / Fnum do
	sol := fsolve(subs(F = x,det) = 0,lambda);
	points := points union {seq}([x,sol[i]],i = 1..5)
end do:
plots[pointplot](points);

The following code, where I have just guessed at some expressions for the functions g[j], works fine:

for j from 0 to 10 do
	f[j]:= a[j]*exp(I*2*j)+b[j];
	g[j]:= exp(I*2*j)+b[j];
end do:
for j from 0 to 10 do
	s[j] := solve({f[j]=0,g[j]=0},{a[j],b[j]})
od;

Notice, however, that I use {}-brackets around a[j] and b[j], not []-brackets, which, if used, produces the following error "Error, (in solve) invalid arguments".

Let A and B be given by

A := Matrix([[x,y],[u,v]]);
B := Matrix([[1,2],[3,4]]);

say. Then,

for i from 1 to LinearAlgebra:-RowDimension(A) do
for j from 1 to LinearAlgebra:-ColumnDimension(A) do
	assign(A[i,j],B[i,j])
end do
end do:
x,y,u,v;
                           1, 2, 3, 4

Is the following of any help?

for i from 1 to nops(L) do
	K[i] := map(D -> LinearAlgebra:-Transpose(D) . L[i] . D,P)
end do:

In Maple 11 one could do something like:

makeList := (p::posint) -> [seq](Matrix(2,2,convert(i,base,p)),i = 0..p^4 - 1):

It is, however, slower than the other proposals.

Please, disregard this post. The previous content was misconceived.

Exporting the plot as a gif-file, opening it in PhotoShop, and there using the eyedropper tool yields the RGB values (179,230,179), which divided by 255 gives (to two decimal places) (0.70,0.90,0.70), consistent with what Georgios Kokovidis finds. But, admitted, resorting to this purely graphical method is a bit crazy when, as Georgios Kokovidis nicely makes one aware of, one can just use showstat(leftbox).

Maybe you know this (if so, just disregard this post), but whenever you see something like x^2 + y^2 + z^2 you should think of spherical coordinates (r,theta,phi), in which the Cartesian coordinates (x,y,z) may be parametrized as

x(r,theta,phi) := r*sin(theta)*cos(phi);
y(r,theta,phi) := r*sin(theta)*sin(phi);
z(r,theta,phi) := r*cos(theta);

They have the pleasant property

> simplify(x(r,theta,phi)^2 + y(r,theta,phi)^2 + z(r,theta,phi)^2);
                                2
                               r 

from which it follows that the surface S of the example you give can be parametrized as

x(theta,phi) := a*sin(theta)*cos(phi);
y(theta,phi) := a*sin(theta)*sin(phi);
z(theta,phi) := a*cos(theta);

For finding a point inside each bounded region, please see my blog post Regions bounded by lines.

Do you mean x^2 + y^2 <= R or x^2 + y^2 <= R^2? The former is what Israel assumes, the latter what Jakubi assumes. I will assume the latter. Because f(x,y,z) = 7z is independent of x and y the integral separates and becomes int(7z,z = 0..m) times the area of the base of the cylinder. So, the result is

> int(7*z,z=0..m)*Pi*R^2;
                           7  2     2
                           - m  Pi R 
                           2         

PS: Reading your post again I can see that a solution is not actually what you ask for. Anyway, here it is.

At present I have coded the following:

• Generating some coefficients (the factor 5 is present to get some interesting plot):

  n := 8:
  gen := rand(10000) - 5000:
  a := [seq(gen(),i=1..n)];
  b := [seq(gen(),i=1..n)];
  c := [5*seq(gen(),i=1..n)];
  

• Generating the equations for the lines:

  eqs := [seq(a[i]*x + b[i]*y = c[i],i = 1..n)];
  

• Plotting the lines:

  plot(map(solve,eqs,y),x = -10..10);
  

  giving, for instance,

But I have no clear idea yet as to how to obtain the coordinates of some points inside those regions. Maybe the fact (unless I am fundamentally mistaken) that all bounded regions are convex polygons is helpful: most probably there exist some powerful procedure in such a case. Hopefully, some hardcore mathematician, which I am not, can help us here.

Two minor comments:

• As written, your linear equations all describe lines passing through the origin, and so they define only unbounded regions. Instead of A[i]*x+B[i]*y = 0, do you mean A[i]*x+B[i]*y = C[i]?
• Disregarding the degenerate cases, where one or more lines are parallel or identical, there are not (8-1)! = 5040 intersection points between pairs of lines, but only binomial(8,2) = 28.