@Carl Love There is viscous damping between the ball and the container, meaning a resistive force which is proportional to their relative velocities. There is dry friction between the container and the floor, and that's responsible for the jerky motion.

The viscous damping enters the equations of motion through the term c*(v₂−v₁). Taking c=0 is possible but in that case the ball will forever oscillate after the container has come to rest.

As to marketability, I don't know; that's above my pay grade :-)

@Preben Alsholm I had completely missed the event_doublecross option in the doucmentation, and that's exactly what is needed here.  Thank you very much for pointing it out.  Now I am happy :-)

@Carl Love Motivated by your comment, I made a model of something similar to sloshing water in a bucket.  The "water" in the model is just a point of mass which bounces back and forth in a container which in turn can slide on the floor.   I think that calculation deserves a post of its own, so I will post it in a new thread.

@mmcdara I have read your suggestion regarding an alternative approach to treating dry friction and I like it quite a bit.  The result produced through it are quite compatible to those obtained through dsolve/events and have the advantage that they involve straightforward numerical solution of a system of differential equations and therefore are more robust than the delicate approach through events.

I will post a solution of a different problem which I solved by applying the friction model suggested by you.

Thanks for the very helpful input.

@mmcdara I am out on travel now. I will return in two days and carefully read your comments and then reply.  I thank you very much for your extensive feedback.

@mmcdara I have tried smoothing out the friction term but have found the results to be unsatisfactory in the sense that even with quite small eps (in the notation of your code fragment), the solution remains noticeably different from that obtained without smoothing. The discrepancy is most noticeable in the plot of v(t).  The original v(t) has flat parts, indicating a motion with constant velocity (sticking to the belt).  In the smoothed out version, v(t) has nothing resembling flat parts.  This is with eps=1e-3.  Changing eps to 1e-4 leads to a premature halt in the solver as in the unsmoothed case.

@acer Thanks for the feedback.  I rarely think of Explore as an animation tool but I can see that it can be very helpful, especially during code development. As to DocumentTools:-Tabulate, I will have to read the docs because I know nothing about it.  

@Carl Love The automatic root finder is a good idea but we should be aware that it can skip over some roots.  For instance, the roots 245.3396184, 481.5561050, 530.0902992, 578.7535346 are missing in what you have shown.

@mehran rajabi The usual way for searching for roots is to plot the function in order to get a rough idea about where the roots are.  Then call fsolve with appropriate ranges to pinpoint them.  In your case we do:

plot(eq, x = 50 .. 200, discont = true, view = -100 .. 100)
fsolve(eq, x = 60 .. 70);
                          68.93033112
fsolve(eq, x = 70 .. 120);
                          106.9740567
fsolve(eq, x = 120 .. 170);
                          151.0624355
fsolve(eq, x = 170 .. 210);
                          198.3826229

@Adjal yassine I am quite familiar with the finite element method and its uses for solving differential equation, especially those that relate to structures. I am puzzled, however, by your statement that you are interested in the stiffness matrix, but not in the solution.  What is the use of a stiffness matrix other than for computing a solution?

Perhaps you have good reasons for writing your programs but those reasons do not come through in your posts.  If you wish to attract people's attentions, you should be more explicit in explaining your goals.

@acer That's excellent.  I have your setsize in my ~/.mapleinit now.  It works even with animations. Thanks!

The presentation in your worksheet makes things more difficult than is necessary.  Consider the following.

1. You have defined the variables y, B, R, R2, r as functions.  To make life easier, I suggest making all of them into expressions.  For instance, replace

y := (x,k,xi,B) -> whatever

with

y := whatever

2. In the definitions of R and R2 replace B(k) with B.

3.  Is the y that enters the definition of r the same y that is defined on the first line in the worksheet?  I don't think so.  And if that is indeed a different y, then you should name it something different.

4. In order to plot any of these functions, you need to specify the value of k.

@CEFOG Look at y = a*x + b.  Can you solve for y/a so as to get y/a = (something that has no a in it)?  Do it by hand, not Maple. If you show me how you do that by hand, I will show you how to do it in Maple.

@Joa The symbolic solution that you have obtained is too unwieldy to be useful.  At least I cannot make it do what you want.  But the problem that you wish to solve is a well-explored subject and much has been written about it.

Specifically, we have a system of differential equations which involve a certain number of parameters.  We wish to determine the parameters so that the solutions fit a prescribed set of data.  Doing that is not very difficult but requires some familiarity with common techniques in numerical computing.  If you have a serious need (such as graduate research) for solving the problem that you have stated, I suggest that you have a look at Notes on Adjoint Methods.  Section 6 addresses exactly what you need.

@Earl Your calculations are correct, however I would organize the worksheet somewhat differently.  You may continue doing it your way, or do it as shown below, whichever makes better sense to you.

Download Slinky2-alt.mw