@macmp91 You can't effectively make an assignment to both B and B[b] in the same code. One assignment will overwrite the other. If you want to have B subscripted by b for stylistic reasons (i.e., because you like the way that it appears in the prettyprinted code), then use B__b (double underscore).

@Preben Alsholm I'm a bit skeptical about the sum of the squares of the residuals (i.e., the vertical distance between point and curve) being the objective function. Kitonum's function is clearly a better fit of the amplitude of the data. DirectSeach:-DataFit comes with 9 stock objective functions, including sum of squares of residuals. I haven't tried any of these yet.

Just eyeballing your graph, it looks like A*x^2*sin(B*x) where A and B are positive constants to be determined from the data. A is approximately 0.015 and B is approximately 280.

@acer I find it curious that all three solutions presented so far converge at a midpoint deciding that the desired value is a root of cos(16*x). This was not a goal of my solution, which I produced independently of the other solutions.

@acer What does _S000002, which is the argument to arccos in two of the cache entries, represent?

Your error isn't "Newton iteration is not converging"; it's "Initial Newton iteration is not converging." I changed the title of this Question to refect this. It may seem like a trivial distinction until you know that both are valid error messages from the numeric solution of BVPs but the steps required to overcome them can be vastly different. In particular the "continuation parameter" technique, discussed below, can only be used to overcome the "initial" error.

@Ronan The commands that both Kitonum and I gave aren't meant to be executed in isolation. They are meant to be executed only after ra has been defined in terms of x by the code given in the original Question.

@tomleslie From my very limited understanding of Maple's numeric BVP algorithm, I would guess that there would never even be a "Newton iteration" if the equations were linear because the problem reduces to solving a system of linear algebraic equations. Thus I try to eliminate the nonlinear terms if that doesn't make too big of a change to the overall system. If my guess is correct, then this guarantees at least that the "kickstart" at CONT = 0 will work. Then this solution provides initial guesses for the Newton's-method algorithm that's applied to the system of nonlinear algebraic equations at the next value of CONT. 

You are lucky that the easy fixes provided by Tom and Preben work. In my experience, it usually takes more (and more subtle) effort with more guesswork to overcome this error.

I just mention this for the benefit of anyone who comes across this thread in the future and doesn't understand why these fixes aren't working for them.

@Rouben Rostamian You're probably familiar with the concept that Maple's infix operators such as + and * have prefix forms such as `+` and `*`.  (Let me know if you're not familiar with that because I'd be happy to explain it.) These are sometimes required, for example, when you want to use the operator with map or if you want to overload it. In the Maple expression A[x] (regardless of the type of A), there is an implied infix operator between the A and the [x]. This operator has rules of precedence and associativity just like all the explicitly visible infix operators. The prefix form of this operator is `?[]`. This operator's action is pretty much*(footnote) as if it were defined

`?[]`:= (A, idx::list)-> A[idx[]]

except that it's builtin.

This operator is documented at ?use and in Chapter 9 of the Programming Guide. This operator can be overloaded so that you can prescribe an action for expressions A[x] for A of types for which no action is prescribed by default. Or you can override the default action. Or you can simply use `?[]` as a prefix form for map, which is what I did in the Answer above.

Footnote: I say "pretty much" because there's some extra subtlety when A is a "naked" expression sequence or when A[x] appears on the left side of :=. These subtleties are documented, albeit rather tersely, at ?use.

@Vortex Oops, I forgot to check for HFloat(undefined) instead of Float(undefined). The H is for Hardware. So, if I check the matrix with

numboccur(M, HFloat(undefined))

then numerous occurences are reported, concurring with the plot attached to your most-recent Reply. The NaN from the "hardware double precision" arithmetic following the IEEE-754 standard is HFloat(undefined).  The NaN from Maple's internally processed "software floats" is Float(undefined). The automatic conversions between software and hardware floats is a subtle issue probably not worth going into right now.

@kencom1 Forget about the journal for the moment. Focus on your code and what I said above. You assigned an expression to psi, so you can't independently vary psi over a range. You did not assign any expression to y. In particular, you did not assign to y an expression containing x, so the commands subs(x= ..., y)  make no sense.

@Vortex Your correction of my code is correct. What I meant was

(i,j)-> R0f(1 + (i-1)*(50-1)/(50-1), 0 + (j-1)*(5-0)/(1000-1))

That first argument is of course identical to i. I write it like this to emphasize that you are choosing 50 evenly spaced points in the real interval 1..50. There's no need for a to be integer or for the width of the interval to have any correspondence to the number of points chosen. So, a and x are given the same treatment.

Btw, I checked for points where there was no solution for y in 0..1, and I didn't find any. Those points, if any, can be counted with

numboccur(M, Float(undefined));

Maple's name for NaN is Float(undefined).

@Scot Gould A physical justification for gravity proportional to 1/r is if the bodies are long thin cylinders perpendicular to the plane of motion. The bodies would need to have length significantly greater than the distances between the bodies.

Having a procedure that computes a matrix's entries can be better than having the matrix itself. What do you want this matrix for?