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.

I cannot understand your equations.  Nor can Maple.  You can help yourself by writing a more readable code,  How many differential equations do you have?  How many unknowns?  I suggest something like:

de1 := diff(f(eta), eta, eta, eta) + ...
de2 := ...
de3 := ...
bc := f(0)=1, D(f)(0)=...

and then

dsol := dsolve([de1, de2, de3, ..., bc], numeric);

@Volker Lehner In your worksheet some of the entries are German text and some others are Maple code.  Therefore I assume that you know how to enter text and code as needed.

You entered my code as text in your document.  But you should have entered it as code. Here I have corrected that and have illustrated how to run the procedure.  My comments are in red.

Berechnungen-1-fixed.mw

@Preben Alsholm That's very nice analysis.  The "magic" happens between ineq1 and inqe2.  The calculation up to and including ineq1 is consistent with y(x)≥ 0. The strict inequality, y(x) > 0, enters in the separation of variables calculation that takes us from iineq1 to ineq2. Quite insightful.

@JohnS I just posted a response to oceanxinlie in which I explained some of the details of the calculations, including how Lcrit was calculated.  The title of the post is "Details".

@oceanxinlie You asked:

Can I know where do you find the elementary function for solving this differential equation, is it very famous?

Your differential equation y'' = a*y^q has quite a simple form, therefore is not difficult to guess that it must have a solution of the form y(x) = A*x^B.

To determine A and B, we substitute that guess into the equation and equate the exponents of x on both sides and arrive at B = 2/(1-q).  Then we equate the coefficients of x on both sides, whence we arrive at the value of A that you see in the worksheet.

The solution y(x)=A*x^B has the property that y(0)=0 and y'(0)=0.  This does not solve your boundary value problem because you want y(0)=0 and y'(L)=b.  But that's alright.  You want L=2945 but let's ignore that for the moment.  Let's find a different value of L, let's call it Lcrit, so that y'(Lcrit)=b.  Since y(x)=A*x^B, this says A*B*Lcrit^(B-1) = b.  We solve this equation for Lcrit and arrive at the value of Lcrit given in the worksheet.

This would be the end of our calculations if your L were exactly equal to Lcrit.  But it isn't.  Your L is quite a bit larger than Lcrit.  But that's not a problem.  It is not hard to see that if a function y(x) solves the differential equation, then the shifted function y(x-k) also solves the equation for any k.  So we pick k = L − Lcrit so that the shifted function satisfies your desired boundary condition when x = L.

We fill the gap on the x axis between the origin and the shifted function by zero.  This is a well-known idea which Preben immediately recognized with his comment about y(x0) and y'(x0).

You also asked:

Can you show me the code of numerical solution?

That's exactly what Preben did in his first example.  He took L=450 which is less than Lcrit, and found the solution numerically.

The usual numerical methods don't work when L is larger that Lcrit, but that's not a problem because in that case we have an exact solution as shown in the worksheet.

In the rest of your comments you say that you had hoped for a solution which is not flat zero on a part of the interval.  As far as I can see, if L > Lcrit, that zero part cannot be avoided.  Mathematically there is nothing wrong with that.  If you are applying the result to some physical model, you will have to see what that is saying about the model.

@mohammad1369 If I understand your question correctly, you want to write A_x instead of A.  You enter A_x in Maple by typing A__x, that is, A followed by two underscore characters and an x.  For example:

A__x := <a, b; c, d>;

@Les Sure, there is nothing preventing from putting that construction inside a procedure.

@dpaddy Before you get busy coding your function in C, see if this can help you:

h := t -> exp(I*t);
plot(Re(h(t)), t=0..2*Pi);

@Earl If the disk's initial angular velocity is small, then it will rock back-and-forth, rather than take off and roll forward forever.  Note that unlike your drawing, I am taking the blue disk to be at the lowest position at time zero.

The transition from rocking to rolling motion is obtained by giving a sufficiently large initial angular velocity that enables the blue disk to rise from the bottom position to the top position.  At the borderline case, the blue disk will creep to the top position and stop there.  (By the way, in your post where you write Pi/2 it should be Pi.)  That borderline case is determined by equating the initial kinetic energy to the change in potential energy.  The expressions for the two energies are available in the worksheet.  Equating the two and solving for omega^2 we get

omega^2 = 8*g*m*d/((2*R^2-4*R*d+2*d^2+r^2)*m+3*M*R^2), that is:

Plugging in the parameters used in the worksheet we get omega^2 = 0.4232804233, and therefore omega = 0.6506000486, which is what you have observed.

As to g=1, that depends on the choice of units.  The value of g may be 9.8, or 980, or 32.2,, or anything else, depending on the choices of the units of mass, length, and time.  Since we are not saying what those units are, there is no reason to be so particular about the value of g.

Regarding your and JohnS's question about what happens after the vertical contact force hits zero, I expect the disk to become airborne and spin freely for a while, and then fall back and regain contact.  I haven't done the calculations for that.  I may look at it this weekend if I find the time.

@max125 The command int~ is  an instance of Maple's element-wise operation.  Look up the help page on "Element-wise Operators".  In short, int~ applies the integration operator to the individual components of a vector.

I must add that vv's answer is the best solution to your question.  That, too, uses int~ to integrate the vector.

Your solution looks good to me.  Here's a slight variation:

restart;
a := < 4*t, sin(t), cos(2*t) >;
v0 := <1,0,0>;  r0 := <0,1,0>;
# calculate v:
int~(a, t) + C;
eval(%, t=0) = v0;
simplify(isolate(%, C));
v := eval(%%%, %);
# calculate r:
int~(v, t) + C;
eval(%, t=0) = r0;
simplify(isolate(%, C));
r := eval(%%%, %);

There is nothing wrong, mathematically, with your equations.  I assume you are entering them incorrectly.  Post your worksheet so that it can be examined.

Click on the big fat green arrow in the window where you compose your message in order to upload the worksheet.