@alaindika5 Options is on the Tools menu.

@Mariusz Iwaniuk For this to work, it requires that the inner integrals done before the improper one do NOT cause the integrand to become explicitly 0. One may not have that knowledge a priori in practical cases.

Also, it requires that the limits of integration on the inner integral that gets moved outward be constant with respect to the variables of integration that become inwards to it.

@Carl Love The bug only seems to affect multiple integrals expressed with a single int. So, a much better workaround than using limit is to use an old-style nested integral:

int(int(int(
    exp(-2*r)*cos(theta)^3*r^2*sin(theta), 
    phi= 0..2*Pi), theta= 0..Pi), r= 0..infinity
)

or to specify the nesting implicitly with foldl:

foldl(
    int, exp(-2*r)*cos(theta)^3*r^2*sin(theta), 
    phi= 0..2*Pi, theta= 0..Pi, r= 0..infinity
)

@brian bovril This way of using Iterator:-Permute (to get the same 2592 permutations) is not significantly faster than the above, but it's simpler and perhaps easier to understand. In particular, it doesn't use the modular arithmetic.

T:= [[A,X],[A,Y],[A,Z], [B,X],[B,Y],[B,Z], [C,X],[C,Y],[C,Z]]:
S:= {op~(T)[]};
                    S := {A, B, C, X, Y, Z}

P:= output((J:= Iterator:-Permute((``@op)~(T)))):
ps:= [seq](ArrayTools:-Alias(P, k, [3]), k= [0, 3, 6]):
Perms:= eval(
    [in J do if andmap(p-> {seq}(op~(p)) = S, ps) then [seq](P) fi od],
    ``= `[]`
):

A bug in Permute causes it to give incorrect results if the elements being permuted are themselves lists. The code above works around that by temporarily converting sublists of the form [a,x] to ``(a,x).

This produces the same 2592 permutations as Kitonum's latest, but done with Iterator.

restart:
T:= Array(0..8, [[A, X], [A, Y], [A, Z], [B, X], [B, Y], [B, Z], [C, X], [C, Y], [C, Z]]):
z3:= {0,1,2}:
Perms:= [
    for p in Iterator:-Permute([$0..8]) do 
        q:= iquo~(p,3); r:= irem~(p,3);
        for k by 3 to 7 while {seq}(q[k..k+2]) = z3 and {seq}(r[k..k+2]) = z3 do
        od;
        if k=10 then [seq](T[[seq](p)]) fi
    od
]:

@BigGeorge Here's a formula for the surface area of an ellipsoid:

#Surface area of ellipsoid with semi-axes a, b, c. 
#Assumes that a, b, c > 0, c < a, c <= b.
SA__Ellipsoid:= (a, b, c)->
local A:= a/c, B:= b/c, s:= sqrt(A^2-1), z:= s/A, k:= A/B/s*sqrt(B^2-1);
    2*Pi*c^2*(1+B*(EllipticF(z,k)/s + s*EllipticE(z,k)))
:

Using this in your case (I divide by 2 because your integral only covers half the ellipsoid):

(a, b, h):= (213/2, 144/2, 46):
SA__Ellipsoid(a, b, h)/2;
        / 
2116 Pi |1
        \ 

     144           /11     (1/2)   71      (1/2)    (1/2)\ 
   + ---- EllipticF|--- 305     , ----- 305      767     | 
     3355          \213           40260                  / 

     (1/2)
  305     

     99     (1/2)          /11     (1/2)   71      (1/2)    (1/2)\
   + --- 305      EllipticE|--- 305     , ----- 305      767     |
     529                   \213           40260                  /

  \
  |
  /

evalf(%);
                          34384.48452

which agrees with the previous result.

I hope that someone can derive the formula for SA__Ellipsoid in Maple by integration.

@BigGeorge You wrote:

• But it still freezes up

There's a difference between "freezing up" and proceeding with a lengthy computation---perhaps one so lengthy that it won't finish in any reasonable time. When I do int(jf, [u= -Pi..Pi, v= 0..Pi/2]), I see it performing a lengthy computation, as evidenced by the changing "Memory: ..." and "Time: ..." values shown in the bottom right of the window.

Doing a preliminary simplification of the integrand helps, and changing the order of integration helps more, but neither help enough to completely finish the integration:

jf2:= simplify(jf) assuming real;
int(jf2, [v= 0..Pi/2, u= -Pi..Pi]); #Note change in u,v order.

The integral with respect to v is done for me, leaving a messy integral with respect to u. 

Unless you have some reason to think that you "should"^[*1] be able to get an exact result, you should settle for a numeric result, the 34384. that you already got (and I've confirmed that by applying evalf to the remaining integral).

• ...and then deletes all of my random variables.

I don't understand what you mean by "delete". If you enter a does it not return the numeric value of a? Are you using the stop sign on the tool bar to stop the integration? You should.

• Also, when I do jf = sqrt(jac[1]^2+...), does it create the Jacobian as a simple equation/sum, or does it somehow confuse it with a vector column or something?

It's a simple scalar expression, the square root of a sum. (By scalar in this context, I mean "not a vector", as you probably know.)

Edit: 
[*1] I realize that you're trying to compute the surface area of half an ellipsoid. There is a formula for that, although it's complicated. So, you do have a good reason to expect that the integral "should" be done. It would be good if Maple could derive the formula. I'm writing a followup on the formula.

@ssgmckv The command for 2D contour plots is plots:-contourplot.

It is clear that the Answers given so far are sufficient to solve the resistance problem for any finite simple graph and any resistances. Thus, you shouldn't post new Questions that are no more complicated than that.

@mehdibgh It's not possible to express a nonlinear function as multiplication by a matrix. This is a mathematical limitation, not a limitation of Maple.

@lcz You wrote:

• Here I have not studied the corresponding relationship between graph vertex index and adjacency matrix index thoroughly. However,  it seems that Rows and columns of the adjacency matrix follow the order given by Vertices.

Yes, that is always true.
 

@dharr You wrote:

• I'm assuming the regular inverse is more efficiently calculated than the Moore-Penrose pseudoinverse.

When doing an exact computation, such as is being done here, the regular inverse is many times faster than the Moore-Penrose.

@bstuan You exclaimed:

• It's amazing how a seemingly complicated problem can be solved so simply!

This "textbook exercise" was highly contrived precisely to initially seem complicated yet have a simple solution. I'm a bit dubious of the pedagogy of that.

@ssgmckv You can't name a procedure 1/h0. You could name it h0 or `1/h0` (note the back quotes).

Technically, you are allowed to name a procedure D(D(F))(0), but I think it's a bad idea to do so. You could use DDF0 or `D(D(F))(0)`. Anything is a legal name if you put it in back quotes.

@bstuan See the comment that I just added to my Answer explaining how to do a "by hand" solution. And note that

exp(f(x)-x^2-1) = exp(f(x))*exp(-x^2-1) 

@ssgmckv If you want to do 3D contour plots of solutions of BVP system with the axes representing parameters and/or endpoint values (such as D(h)(0)), then see this Post: "Numerically solving BVPs that have many parameters"