@nm Oops, sorry.  I didn't notice the ".1"     :(

@sarra The worksheet that you have started has several issues which makes it difficult to work with.  Consider the loop:
for i from 1 to 4 do
    whatever;
end do;
After that loop ends, the symbol i has the last value assigned to it in the loop, which is 5 in this case.  From that point on, if you enter something like b[i] in your worksheet, it evaluates to b[5] which is generally not what you want.

There are several ways of getting around that.  The easiest way is not to use for-loops at all.  Look at finite-difference.mw to see how to do that.  You should be able to extend it easily to complete the solution of your problem.

Remarks:

• I haven't checked the correctness of your mathematics at all.  Check!
• I have changed beta1, beta2, and gamma to beta[1], beta[2], and beta[3].  The first two are cosmetic.  The third one is important, because gamma is a predefined constant in Maple.  You should avoid using it.
• I suggest that you use Maple document format for your worksheet and 1D notation for the input because that makes life easier for problems such as yours.  Instructions are given here:
  http://userpages.umbc.edu/~rostamia/math481/config/maple.html

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