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

Carl, what I am asking is why is it that g[12] has any meaning at all?  Where is such a construct defined/documented?

@vv Your construction puzzles me.  Let's look at a simpler case:

g := x -> x^2;

Then g(5) returns 25.  I also see that g[12](5) returns 25 but I don't understand why.  What is the meaning of g[12]?  Would you please explain?

@Volker Lehner Short answer: To enter new material inside a procedure, type Shift-Enter.

For anything other than toy programs, I must echo acer's suggestion of using 1D input plus the "Worksheet" instead of the "Document" mode.  You will find instructions https://userpages.umbc.edu/~rostamia/math481/config/maple.html

As far as I know, you cannot convert your existing worksheet.  Configure Maple as instructed, then open a fresh new worksheet.