@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.

@dharr Thanks for your answer and a thumb up!

A curious side note:

Before I posted my question I did exactly what you did, however, in my actual problem I had

eq:=((4*n-1)*x + 4*n + 1)*x^(2*n) = 1 - x;

where the coefficient of every n is twice that in the equation which I posted.

As I was posting the question, I noticed that all coefficients of n are even, therefore I replaced every n by n/2 without trying to run the reduced version through Maple.  Curiously, the calculation with the original version fails:

eq:=((4*n-1)*x + 4*n + 1)*x^(2*n) = 1 - x;
solve(%,x);
asympt(%,n);
Error, (in asympt) unable to compute series

I don't understand why it should work for n but not for 2*n.

@permanoon123 In his very first response, tomleslie pointed out that the second derivative with respect to x is written D[1,1] in Maple.  Your D[2] means the first derivative with respect to t, but that's not what you want.

Additionally...

I see that you corrected the sign of the fourth order derivative term relative to your initial post.  At the same time you also changed the sign of the second order derivative term.  That's bad!  The sign of the second order derivative should be positive.

Are you making up random PDEs?  Don't do that!

@Chrono1 I have added a simplified front-end to the worksheet which does what you are asking for.  Specifically, the original geodesic_points() expects the initial segment PQ of the geodesic to be specified by (a) the face number in which PQ lies; and (b) the edge numbers in that face on which the points P and Q lie.

Since all edges and faces on a regular polyhedron are equivalent, there is no loss of generality in always taking the line segment PQ on Face 1, and the point P on edge 1 of that face.  Furthermore, instead of specifying the point Q, we may specify the slope of the line PQ relative to edge 1 and let Maple find the point Q.

This alternative way of specifying the segment PQ is implemented in the proc geodesic_points_alt() at the end of the attached worksheet.  See if that meets your needs.

Download cube-ver2.mw

@tomleslie Technically there is no justification behind applying a naive finite difference method to solving a PDE with discontinuities.  Anything we get out of that is a "gift" rather than an entitlement.  In this case Maple's finite difference algorithm picks up a reasonable approximation to the exact solution, that is u(x,t)=Heaviside(t - x/4), as we see in the graph below where green and red colors correspond to the exact and approximate solutions, respectively.

I calculated the approximate solution by setting spacestep=0.001 in pdsolve().  A provably good approximation may be obtained through shock-capturing methods, but that's a highly specialized subject.