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.

If you define the pointers from the top level downwards, then everything works as you expect. In your example, b is at the top level, and m and a are both at the second level. Please execute this:

restart:
b:= m[a,1]:
m:= <1,2; 3,4>; #easy 1-D way to define small matrices
                                      [1  2]
                                 m := [    ]
                                      [3  4]

a:= 1:
b;
                                      1

a:= 2:
b;
                                      3

In the implicitplot command, it is not the purpose of the arguments x= ..., y= ... to set the length of the axes. Rather, their purpose is to define the mathematical search space in which solution points (x,y) to the equation are searched for. The view option is used to set the (logical) length of the axes: view= [-10..10, -10..10].

The rank problem is only due to poor "conditioning", i.e., extremely small singular values of the matrix generated from the input data being erroneously interpretted as zeros. (A much-more-detailed explanation of this phenomenon is in Acer's Answer.) This can be avoided by shifting the data toward the origin (and optionally scaling it). In the code below, each data point is converted to a pair of z-scores before doing the PolynomialFit. Then the returned function is scaled and shifted back to the original units. This trick is "safe" in that it won't work if the raw input data is truly of substandard rank.

restart
:
Rok:= Vector([2013, 2014, 2015, 2016, 2017, 2018], datatype= hfloat):
TrzbyCelkemEmco:= Vector(
    [1028155, 1134120, 1004758, 929584, 995716, 1152042]*1e-6,
    datatype= hfloat
):
St:= Statistics:
SM:= St:-PolynomialFit(
    3, 
    (Rok -~ (mu_x:= St:-Mean(Rok)))
        /(sigma_x:= St:-StandardDeviation(Rok)),
    (TrzbyCelkemEmco -~ (mu_y:= St:-Mean(TrzbyCelkemEmco)))
        /(sigma_y:= St:-StandardDeviation(TrzbyCelkemEmco)),
    x,
    summarize,
    output= solutionmodule
):
#Shift back to original positions:
KubickaTrzby:= (
    mu_y+sigma_y
        *subs(x= (x - mu_x)/sigma_x, SM:-Results('leastsquaresfunction'))
):
printf("\n\nKubicka trzby: %a", evalf[7](KubickaTrzby));
printf(
    "\nCoefficient of determination (r-squared): %a",
    evalf[4](SM:-Results('rsquared'))
);
printf(
    "\nR-squared (adjusted): %a", 
    evalf[4](SM:-Results('rsquaredadjusted'))
);

Summary:
----------------
Model: -.62638474-1.8421266*x+.75166169*x^2+1.3323389*x^3
----------------
Coefficients:
              Estimate  Std. Error  t-value  P(>|t|)
Parameter 1   -0.6264    0.3333     -1.8795   0.2009
Parameter 2   -1.8421    0.6664     -2.7643   0.1097
Parameter 3    0.7517    0.3039      2.4731   0.1319
Parameter 4    1.3323    0.4316      3.0871   0.0909
----------------
R-squared: 0.8874, Adjusted R-squared: 0.7185

Kubicka trzby: 171.5772-.8463921e-1*x+.6461131e-1*(.5345225*x-1077.330)^2+.1145251*(.5345225*x-1077.330)^3
Coefficient of determination (r-squared): .8874
R-squared (adjusted): .7185

Here is the plot with several more "bells & whistles" added:

#Include a margin for axes' view option:
Margin:= proc(min, Max, margin)
local S:= margin*(Max-min);
    min-S..Max+S
end proc
:
plot(
   [<Rok | TrzbyCelkemEmco>, KubickaTrzby], 
   x= `..`(((m,M):= (min,max)(Rok))),
   style= [point, line], color= [red, "DarkGreen"], thickness= 3,
   legend= [`raw data`, `\ndegree-3\npolynomial\nfit`], 
   legendstyle= [location= right],
   symbol= soliddiamond, symbolsize= 12,
   axes= boxed, axesfont= [Times,12], size= [1300, 800],
   labels= [
       'Rok', 
       typeset("Trzby Celkem Emco ", `&times;`, 10^`-6`)
   ], 
   labelfont= [Helvetica, Bolditalic, 16],
   labeldirections= [Horizontal, Vertical], 
   view= [Margin(m,M,.1), Margin((min,max)(TrzbyCelkemEmco), .1)],
   title= "Kubicka Trzby\t\n", titlefont= [Times, Bold, 24],
   caption= typeset("\nData fitted to ", evalf[6](KubickaTrzby)),
   captionfont= [Times, 14]
);

It's a package command. So, either use DEtools:-phaseportrait (my preference) or with(DEtools). 

I'd recommend against trying to reassign ||; doing so could easily make other things not work. But it's very easy to define an infix operator &|| like this:

`&||`:= (R::seq(algebraic))-> `if`(nargs=0, infinity, 1/`+`(1/~R)):

And call it like this:

1 &|| 2 &|| 3;
R1 &|| R2;

I wrote it to take an arbitrary number of operands (or arguments) with null input defaulting to a resistance of infinity (which is the identity operand for this operation). Feel free to change that.

It's safe to overload most predefined infix operators if you can come up with a way for it to recognize from the types of the operands that it's supposed to use its newly defined behavior instead of its original behavior. But if you want an infix operator without any predefined meaning, it must start with &.

@suvetha2000 Okay, I found it: It takes just two commands to convert a high-order linear ODE (or system thereof) to a first-order system in matrix form: DEtools:-convertsys and LinearAlgebra:-GenerateMatrix. The first converts it to a system of algebraic equations (and doesn't actually require linearity); the second extracts the coefficient matrix (and additional vector for the nonhomogenous case).
 

Download ODEtoMatrix.mw

There are at least three reasons why you can't or shouldn't use unapply in this situation: 

1. The most important reason is that the variable that you're trying to unapply unfortunately is used as both the independent variable (naturally) and the variable of integration. I think that this is shoddy on Maple's part. In first-year calculus we teach students not to do that! The correct integral command is Intat rather than Int, and dsolve usually does this correctly. Fortunately, there is an easy way (given below) to work around this. The Int(1,0) that you got is nonsense that results from substituting t=0 into the integral, as can be easily seen.
2. The next reason is that you can't use unapply on a sequence, such as the sequence of two solutions that you got from dsolve. At the minimum, you'd need make the sequnece into a list first.
3. (This one is only a minor stylistic consideration.) The third reason is that you have equations rather than expressions. The more-usual thing would be to unapply just the rhs (right-hand side) of the equations.

A way to get around these problems is to give an initial condition to dsolve. It doesn't even need to be that specific. For example, r(0) = r0 will work:

dsol:= dsolve({diff(r(t), t) = r(t)*(1-r(t)^2) + 2*r(t)*cos(t), r(0)= r0});
R:= unapply(rhs(dsol), [t, r0]); #Make it a function of both t and r0.

This symbolic solution doesn't have much practical value for numeric evaluation. You'll need to use evalf to get numeric integration after giving numeric values for t and r0.