@Ramakrishnan Your code does produce the required result but it calls for manually prepending and appending left and right extensions to the given data.

The solution that I proposed has the slight advantage in that it produces the extensions automatically; no user intervention is required.

I don't understand the boundary condition Dt(theta(0, t)) = 0.  I only guess that Dt is defined in one of the packages that you have loaded, and it looks as if it may mean the derivative with respect to t which in the standard Maple notation would be written as D[2](theta)(0,t) = 0.

If that is so, then I am afraid that the problem is not well-posed.  Even in the case of the standard heat equation the boundary condition D[2](theta)(0,t) = 0 looks very odd to me.  Do you really mean it?

@Ramakrishnan The objective is to produce a function which is constant outside of the prescribed range.  Your function continues growing forever.

@9009134 You are still focusing on the contents of the file 'u(y)-y.dat' .  That's the wrong place to look.  Before you put anything in that file, you should verify that there is something useful to put in it.  Did you watch the output of lprint(tt,u) on the screen? Is that what you expect? Probably not.

Fix that error first, and then we will talk about why 'u(y)-y.dat' is empty.

@9009134 Forget about the file 'u(y)-y.dat' for now. Did you do what I suggested?

Temporarily replace the fprintf line with lprint(tt,u) so that you can inspect the values of tt and u.  You will see that u does not evaluate to a number. Track down the error.

The Q[..,3]:=Q[..,3]/Q[1,3] calculation adjusts the third eigenvector.  In general to adjust all columns in a matrix Q so that their first components are 1 (if possible) we may do:
Q := Q . DiagonalMatrix(1/~Row(Q,1));

@leafgreen As an alternative to Kitonum's correction, you may consider defining f as an indexed variable, as in:

f := Array(antisymmetric, 1..3, 1..3,
    (i,j) -> -m[i]*m[j]*(x[i]-x[j])/((x[i]-x[j])^2+(y[i]-y[j])^2) );

Then f[1,2] will produce the x component of the gravitational force.  You need to define another function to produce the y component.

I should point out that your force formula assumes that the gravitational force is proportional to the inverse of the distance between the particles.  Do you really mean that?  Newton says that the gravitational force is proportional to the square of the inverse distance.

In either case, I am afraid that your worksheet takes you in a direction which is not suitable at all for implementing Verlet's algorithms.  What happens later when the particles have moved?  Their distances change, and the forces need to be recomputed.  The current setup is not quite suited to that.

Someone needs to grab you and turn you around by about 90 degrees to point you in the correct direction.  Your professor should be able to help with that.

@Franzs If you want the axes to rotate along with the object, define P1 as before, and then do:

p := Plot(P1, view = [-40 .. 40, -40 .. 40, -40 .. 40],
    color = "Orange", thickness = 2,
    viewpoint=["circleleft", frames=60], axes=normal);

I don't understand what you mean by "without the names of its spatial coordinates".

@rrbaldino  You keep saying that the second answer is the solution although you plainly see that it fails to satisfy your initial condition.

Think about this:  If I give you several candidates for solution, how would you pick the right one?  Do you know a way other than checking that the differential equation and the initial conditions are satisfied?  Is there some other magical way?

Let me repeat what I said before:  Your problem is ill-posed.  The way to fix it, is to shift the time of the application of the delta function to somewhere after the initial time.

What did you get for the critical point?  My calculation gives T = -45926.73299, v = 0.2514345875.

But this smells like a question coming from thermodynamics, and therefore that temperature looks a bit suspicious to me.

@rrbaldino  Did you plot the two "solutions" as I suggested?  Neither of them satisfies the initial condition dy/dt=0 at t=0 that you have asked for.  One or the other answer that you have obtained may be a solution to some problem, but certainly not to your problem.

If anything, Maple should warn that both answers are wrong.

@vv What you have here is a tour de force of recursive programming.  Awesome!

I think that you are asking how to convert ln(x)/ln(10) to log₁₀(x). So far as I can tell that's not possible—Maple automatically converts logarithms from any base to base e.

For displaying purposes you may temporarily prevent the conversion by putting the desired expression in quotations, as in 'log[10](x)'; which yields log₁₀(x). But if you refer to that result later, it will revert to natural logarithms.

@nikkarkian In mathematics, "unknown" is what is to be solved for.  In your equation, V is not "unknown" because you are not solving for it.  The unkown is q(t), so all you need is:

EulerLagrange(L, t, q(t));