@Mariusz Iwaniuk That's excellent. A thumb up!  I wonder whether it's possible to recover the original integral through differentiating the answer.

@tomleslie No, I didn't mean to.  I played around with the equation a bit and, among other things, changed the sign of the exponent, but then I forgot to restore it.  Mea culpa!

That said, the function exp(t*x) grows astronomically large very quickly.  I can't think of any realistic problem leading to such a term.  I like the problem better with the minus in the exponent.

@kambiz1199 A complete answer to your question depends on the details of your worksheet that you have not shown, so people are just guessing what you may have in there.

Here is a complete worksheet where for simplicity I have taken the N matrix to be 2x4, and have taken some random expressions for N[1] and N[2].  If this works for you, then you should be able to extend it as needed.

Download mw2.mw

@arshl  We apply dsolve to find the solution v[3](t) of the differential equation.  The solution depends on many parameters such as p1, R, k11, a11, ro, cp2, etc.  There is nothing special about p1.  Maple is quite happy solving differential equations involving parameters.  Here is a simple example:

de := diff(x(t),t,t) + omega^2 * x(t) = 0;
dsolve(de, x(t));

The solution x is a function of t, so we write it as x(t).  It also depends on the parameter omega, but we don't write x(t,omega).

@tomleslie Thanks for pointing this out.  I should have noticed that the integral symbol returned by dsolve is colored gray, indicating an inert Int.

@arshl In my previous reply I said that all "pi" should be "Pi".  You haven't made that change yet.  You should.

Looking further down your code, I don't understand the meaning of 
m1(p1, my(t)) := (1/2)*b11*ro*pi*R^2*(my(t)*cp2/(R*p1)-cp3)/p1;
What are your attempting to do there?

Perhaps it would be easier if you explain just the mathematics of your problem (no Maple).  Then someone may suggest ways to convert that mathematics into Maple.

Your say you are applying dsolve, but you haven't shown the code that produces the result that you have shown.  Showing only a part of your worksheet is not very informative.  Post your entire worksheet as an attachment.  How to do that? In the window in which you edit your message before posting, note the big green arrow near the top edge.  Click on that arrow to attach your worksheet.

But...

Before you do that, change all occurrences of "pi" in your worksheet to "Pi".  They are not the same thing.

@escorpsy

Download mw.mw

@mmcdara I pointed out the ambiguity of the definitions of c_1 and c_2, but farazhedayati said that he has reasons for treating xbar as an independent variable, so I just implemented what he is asking.  It's up to him to make sense of the result.

@farazhedayati To add the c_2 values, we do

add(c_2(my_x, my_xbar, i), i=1..n);

and we get 1.479074043.

The definition
    c1:=diff(b, xbar);
is ambiguous.  To illustrate that in a very simple case, consider the function f = m + 2*x + 2*y, where m is the mean of x and y.  You ask, what is the derivative of f with respect to m?

Answer #1:  Clearly df/dm = 1.

Answer #2:  By the definition of m we have m = (x+y)/2, therefore x+y=2*m, that is, 2*x+2*y=4*m.  Consequently, f=m+4*m=5*m, and therefore df/dm=5.

Which answer would you accept?

This alternative version provides a better demo.

Download stretch-octahedron2.mw

@tomleslie What I was looking for was an asymptotic estimate like the one in dharr's answer.  Nevertheless, thanks for your input and a thumb up for that!

@vv Thanks for your observation.  I don't quite see how you did that by hand.  I will have to look at it more closely.

Added later:

Aha! Got it! Here is how it goes.  We are looking at the equation
  [(4n-1)x + (4n+1)]x^{2n} = 1 - x.
Letting z=1-x the equation takes the form
  [8n - (4n-1)z] (1-z)^{2n} = z.
As n grows, x goes to 1, and therefore z goes to 0.
Consequently, the second term in the square brackets
is smaller than the first, so we drop it and arrive at
   8n (1-z)^{2n} = z.

For small values of z we have 1-z =  exp(-z), therefore
  8n exp(-2nz) = z,
which we rearrange into
  2n z exp(2nz) = 16n^2.
Letting 2nz = y this becomes  y exp(y) = 16n^2, whence from
the definition of LambertW we get y = LambertW(16n^2).
Therefore z = 1/(2n) LambertW(16n^2), and finally
  x = 1 -  1/(2n) LambertW(16n^2).
 

@acer Thanks for your comments, and submitting a bug report.