@tsunamiBTP 

Maple wouldn't yield the result = 1.862645149e-08 because the correct result is 1.931384967026*e-8 as I showed before.

As to why Maple produces 0. by default, consider this.  We have K = 25888158.4218766, therefore within 10 digits of accuracy, we have K^2 = 6.701967464e14.  Therefore within 10 digits of accuracy K^2 - 1 is the same as K^2 because taking away 1 affects digits beyond the level of accuracy that are being carried.  Consequently we get K = sqrt(K^2-1) = 0. which is the correct answer if we perform floating point calculations with 10 digits of accuracy.

If you increase the digits of accuracy, then you get the correct result of 1.931384967026*e-8, as I showed before.

Aside:

The problem that we are faced in that calculation is due to subtracting two large and almost equal numbers.  The usual way to avoid that is to note that

K - sqrt(K^2-1) = (K - sqrt(K^2-1)) * (K + sqrt(K^2-1)) / (K + sqrt(K^2-1))

which simplifies to 1 / (K + sqrt(K^2-1)).  This expression is algebraically equivalent to the original, however it now involves adding, rather than subtracting, two large numbers, and therefore it avoids the earlier numerical truncation issue.  With Maple's default 10 digits of accuracy we see that:

restart;
K := 25888158.4218766;
1 / (K + sqrt(K^2-1));
    1.931384967*10^(-8)

Many years ago George Pólya wrote a book titled How to Solve It. In it he says:

If you cannot solve the proposed problem, try to solve first some related problem. Could you imagine a more accessible related problem?

Which is the simplest finite difference problem that you know how to solve? Try that one first before you attempt this.

@vv By my note I meant to demonstrate a curious construction rather than a serious solution.  Your solution is definitely the right answer to the question.

@ProfJP I am guessing that what you are asking is for the values of x(tt), y(t), x'(tt), y'(t), x''(tt), y''(t), evaluated at a prescribed sequence of times, and where the results are stored as rows of a matrix A.

If that is correct, then start with dsol as I explained earlier, and then do:
dat := [seq(eval([
    x(t), y(t), diff(x(t),t), diff(y(t),t), rhs(de1), rhs(de2)
], dsol(s)), s=0..1, 0.1)];
A := Matrix(dat);

Here time goes from 0 to 1 in steps of 0.1.  Change as needed.

Alternatively, tomleslie's solution may be used for that but I am unable to tell where your difficulty lies from reading your message.  Try to be more clear in what you write.

Cramer's formula shows that each unknown in the solution of a linear system is the ratio of two determinants.  The determinant of an n×n symbolic matrix is the sum of n! terms.  But since
    170! = 7257415615307998967396728211129263114716991681296451376543577798900561843401706157852350749242617459511490991237838520776666022565442753025328900773207510902400430280058295603966612599658257104398558294257568966313439612262571094946806711205568880457193340212661452800000000000000000000000000000000000000000
you will have to wait a bit for your calculation to end.

Preben, that is a very nice construction.  It is particularly interesting that you are able to fit two boundary conditions to a first order ODE.

A slight variant with neater algebra is provided by the ODE diff(y(x), x))^2 = 1-y(x)^2, that is

whose general solution is y(x) = sin(x+c), as well as y(x)=±1.  By piecing together these functions we may fit a solution to the boundary conditions y(0)=a,  y(L)=b (where -1 ≤ a, b ≤ 1).   To guarantee the construction, we need L ≥ 2π.   The solution obtained this way is not unique in general.

For illustration, let's take y(0)=1/2,  y(3π)=sqrt(3/2).  Then it turns out that the function

is a solution for any value of the parameter t in the range 0 to 5π/6.   This animation shows what these solutions look like for varying t:

bvp-for-first-order-ode.mw

@Earl I suspect that the numerical difficulty in solving these equations is due to expressing the solution curve as a function of x.  We should be able to express the solution curve parametrically, as in <x(s), y(s), z(s)> where 0 ≤ s ≤ 1, together with the constraint x'(s)^2 + y'(s)^2 + z'(s)^2 = constant.  I spent a couple of hours today trying to figure out how to do that but did not get anywhere.

If I make any progress in that direction, I will post again.

@Carl Love Thanks for this.  I had seen this trick but had completely forgotten it.  I just uploaded a modified worksheet which includes a section that computes the travel time following your suggestion.

@Yee Voon The letter "I" represents the imaginary unit, that is, sqrt(-1), in Maple.  Use some other letter in your calculations.

@lucaud Your question refers to "inverted pendulum" without saying what it is that you want to do with it.  Because of that we have at least three different proposed ideas:

• My model of stabilizing an inverted pendulum through vertical oscillations of its pivot;
• Thomas Richard's suggestion of stabilizing an inverted pendulum attached to a horizontally moving cart;
• The photo that you posted in reply to Thomas Richard's suggestion which shows something totally different.

To get a useful response you should formulate your question much more clearly.

@MrMarc 

As I said earlier, your system is guaranteed to have the solution {x=0, y=0, z=0}.  If you want a nonzero solution, then you will have to impose the extra condition -100*b*p-p*r+100*a+100*p-r = 0 on the coefficients.  If we solve that condition for r, we get r = (100*(-b*p+a+p))/(p+1).

Plugging in the numbers that you have given for a, b, and p we get r = 13.28571429 which is close to your value of r = 13.29.  That's why it appears to you that you are getting a well-defined solution but that's only an illusion.

You should know that if r is anything other than what I have shown, then the only legitimate solution is {x=0, y=0, z=0}.  If r is exactly what I have shown, then x can be anything, while y = p*x.  For example, with your choice of p = 2.5 and x = 215054, we get y =  p*x = 537635 as you have noted.  But you could taken any other x, let's say x = 1234, in which case you will have y = 1234*p = 3085, and that would have been just as good a solution.

@Earl The change of variable from alpha to t in vv's calculations is not essential.  It's possible to do the calculation directly with alpha.
 

Download projective-alt.mw

@vv Thanks for the pointer to showstat(log) which clearly shows the use of the optional bracketed arguments in Maple's log[b](a) procedure.

What you have shown works for me.  Try re-executing your worksheet.

Aside: Where you have said PLOT(...) I would have said display(...).

@vv I looked up several help pages pertaining to Maple's procs but was unable to find a reference to the f[u](x) construct.  It must be hidden in an obscure place.