@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

@Earl  Page 212 of that book reduces the solution of the problem to solving a system of differential equations for the unknown functions r(m*) and F(m*)  [ equations (3) and (4) ]  subject to the boundary conditions in (2).

The rest of that page muddles things up by attempting to establish an analogy with harmonic oscillations.  The analogy is quite confusing since there is no time-dependence in this problem—we are looking at the steady-state rotation of a stretched spring (unvarying length) about an axis.  There is not much to animate here.

Here I will show how to solve the equations noted above in Maple.  The note at the top of page 213 does the same thing, therefore the confusion of page 212 was not necessary at all.
 

Download mw.mw

@imparter I know that alpha is a parameter that you want to vary.  The problem is, the code that you have shown has a few dozen of ohter symbols whose meanings are known only to you and the author of the book that you are reading.  I can only guess what those symbols mean.  In fact, I can only guess what the code is attempting to solve because you have not even provided me with the mathematical statement of the problem.

That said, I have made a best effort to guess what you have in mind and have produced this worksheet:  fixed-alpha.m which plots the solution curve at t=0.2 as you have asked, but for alpha=1 only.

To plot the curve for other values of alpha you have two choices:

Choice #1. Copy and paste the entire code (but not the "restart" line) in the worksheet (thus doubling the number of lines in the worksheet).  Set the value of alpha in the second half, execute the entire worksheet.  Then you will have two graphs, the original one (with alpha=1) and the second one with whatever value of alpha you want. 

Copy/paste some more to produce additional graphs with other choices of alpha.  Then you may produce a composite graph by applyting plots:-display.

Choice #2. Instead of copying/pasting, you may wrap the entire code in a proc which takes alpha as a parameter.  Then you may execute that proc any number of times with the desired values of alpha.  This second approach is the right way of doing it but it needs some elementary knowledge of Maple programming.  I cannot tell how much of Maple programming you know.  If you don't, then go with Choice #1.

@imparter The plot the curves based on the code that you have shown requires reading the book to understand what the various variables mean. But I don't have the book. Sorry.