Just include fclose("input.txt") in the procedure.

I think that in most languages, most opened files stay open until they are explicitly closed.

Here is a much-faster and slightly more accurate way to solve the BVP by finite-difference methods. This is a method that I just came up with: Start by approximating the function value at the midpoint of the interval. So now we have two subintervals of length (c-a)/2. On each iteration:

1. Divide all subintervals in half.
2. Approximate the function value at each odd-numbered node (the s[k]) using the even-numbered nodes (s[k-1] and s[k+1]) adjacent to it, which have already been approximated. The odd-numbered nodes are the midpoints of the new subintervals.
3. Use fsolve on the entire set of finite-difference equations to refine the approximations (both odd and even nodes), using the existing approximations as initial guesses to fsolve.

This method is limited by the fickleness of fsolve. In actual practice, I'd code a Newton's method solver to take the place of fsolve for this. Nonethless, the worksheet below gets a 127-point solution in a few seconds with a maximum absolute error of 6e-7, (or 9 significant digits). That's pretty good!
 

Download BVP_FDM.mw

This is how I like to do a finite product (and it works just as well for add and seq as for mul) that uses all indices except for a specified one. The specified one in this example is i.

i:= 3:  n:= 5:
mul(1/(x[i]-x[j]), j= [$1..i-1, $i+1..n]);

This works even when i is 1 or n.

Another version:

i:= 3:  n:= 5:
mul(1/(x[i]-x[j]), j= {$1..n} minus {i});
 

Your ode is a quadratic expression in diff(sr(r), r$2)---it's highest-order derivative---and thus it's really two odes, one for each solution of the quadratic. Once you isolate these two branches, it's trivial and very fast to solve each with Maple's BVP solver available through dsolve(..., numeric).

Your solution is sufficiently accurate (shown below), although your technique is doomed to be slow due to symbolic expression swell.
 

Download 2branchBVP.mw

Just do

plots:-odeplot(F, [x(t),y(t),z(t)], t= 0..2);

You can change the range 0..2 to any finite real interval.

We've had an interesting discussion here exploring the purely mathematical and computational details of fitting a polynomial to historical data. But the practical details are that you're trying to model annual sales of oatmeal products. This polynomial model (or any polynomial model) is utterly worthless for predicting future sales: They cannot model the cyclic variation that's evident in your data. Assuming that the cycle is not just randomness (and there's really not enough data to assume even that), it appears to have a period of four years and no growth trend. A sinusoidal model would be more appropriate, but I don't think there's enough data to fit it.

Taking a wild guess: Oatmeal sales are higher in years with Winter Olympics when your company sponsors "home-country hero" athletes who eat a lot of healthy oatmeal. So, find some effective way to advertise in the other years.

Your proposed solution can be immediately seen to be wrong for two obvious reasons:

1. W is a function of two variables, so W' is meaningless.
2. Your solution doesn't contain the derivative of Y.

D[1](w) refers to the first partial derivative of w with respect to w's first parameter. D[2](w) is the first partial derivative with respect to w's second parameter.

Use

subs(s1= s, s2)

One way would be to "zip" it. Let me know if you need more details.

We can vary the complex argument (the t in exp(I*t)) and plot an arbitrary number of roots color coded by that argument.

p:= add(seq((lambda/3)^k/k!, k= 0..4) *~ [1$3, -1/2, 1]);
pi:= evalf(2*Pi):  step:= .01:
plot(
    [seq([Re,Im]~([fsolve(p=exp(I*t), complex)]), t= 0..pi, step)],
    style= point, color= [seq(COLOR(HUE, t/pi), t= 0..pi, step)],
    symbolsize= 4,
    labels= [Re,Im], axes= boxed, size= [900, 700],
    scaling= constrained,
    caption= typeset(
        "\n\t\t", 4*trunc(pi/step), " complex roots of ", abs(p)=1, 
       "\n\tcolor coded by the argument of the polynomial"
    ),
    captionfont= [Times, 14]
);

This gives a clue about those (apparent) circles in the corners: Those roots are uniformly distributed among the input arguments.

Another fascinating property is that the plot is not symmetric with respect to conjugation when you consider the color.

I think that you've completely misinterpreted the purpose of the optional third argument to BinarySearch. The only case where it's useful is where it's the same ordering operator with respect to which the list is already sorted. In a Comment that I wrote to a recent Post, I posted a procedure that works identically to ListTools:-BinarySearch. In the usage examples below that procedure, I showed examples of how those optional arguments could be used. Here's the Comment https://www.mapleprimes.com/posts/212028-Binary-Search-Algorithm#comment203958 

Hints:

1. The price of one additional hour of leisure is the same as the wage that you'd've received had you worked that hour instead. So a decrease in wage is a decrease in the price of leisure.

2. Choosing not to work, i.e., choosing leisure, decreases the labor supply.

3. Your Question can be very easily researched with Google searches. In particular, relevant graphs can be easily found with a Google Image search.

If you find some graphs that you'd like to reproduce in Maple, let me know.

I don't want to do your homework for you, so perhaps this closely related generalization will help you:

limit(c*int((a*sin(x)+b*cos(x))*exp(c*x), x= 0..infinity), c= infinity)    
    assuming c>0; 

The following worksheet is very likely too sophisticated and clever for your direct usage. That's good, because I don't want to do your homework for you. But hopefully you'll learn things from this that you can use. These procedures can be used with exact arithmetic, decimal arithmetic, or even symbolically (i.e., with nonnumeric or partially numeric input).

This should work in any Maple version from the last 10 years or so. If you have any trouble or any questions about the code, let me know.

Download QuadratureRules.mw

If you only want the dependent function values at specific values of the independent variable, then specify output= Array([ list of values ]) in the dsolve command. For your case, that'd be output= Array([$0..10]). If you do this, then the range option serves no purpose and should be omitted,

Manually adjusting the step lengths is generally not an effective strategy if you want high-accuracy results. But using output= Array will help dsolve to choose the step lengths most efficiently to get just the results that you want. Note that all of the "modern" IVP methods used by dsolve (as opposed to the classical subset of methods) use variable step lengths.

If on the other hand setting the step length is more important to you than accuracy (likely for some pedagogical reason), then use a classical method, such as method= classical[rk4], set the stepsize (e.g., stepsize= 1), and specify output= Array(...) as above.