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.

DEtools:-phaseportrait is restricted to rectangular coordinates. You can get around that by using plottools:-transform to explicitly transform the cordinates to polar after it is passed the plot structure from phaseportrait. This causes some distortion of the arrow shapes, but it's not bad. Finally, the plot structure is wrapped in plots:-display so that we can get polar axes.

If you send the trajectories around the circle again by using t= 0..4*Pi (or more), you'll see that the magenta curve (the outermost) is a limit cycle for all of them---like a lopsided cardioid. But the arrows for the 2nd and subsequent revolutions are not the same, so the plot gets cluttered.

In the code, I use t instead of theta simply to reduce the length of the code.

plots:-display( #needed to get polar axes
    plottools:-transform((t,r)-> r*~[cos,sin](t))( #transform to polar
        DEtools:-phaseportrait(
            [diff(r(t), t) = r(t)*(1-r(t)^2) + 2*r(t)*cos(t)], 
            [r(t)], t= 0..2*Pi, 
            `[]`~(r(0)=~.25*~[$1..4]), #inits
            linecolor= [blue, DarkGreen, black, magenta],
            numpoints= 1000
        )
    ), 
    axiscoordinates= polar, size= [1000$2], thickness= 0
);

Here are two more ways:

fsolve(Statistics:-CDF(Normal(mu,10), 500) = 0.05, mu= 500..infinity);
                          516.4485363

fsolve(Statistics:-Quantile(Normal(mu,10), .05) = 500);
                          516.4485363

In the worksheet below, I show how to get the symbolic solution for part d (c = 2.5) by eigenvalue methods. You'd still need to fill numerous of the step-by-step details.
 

Download MatODEs.mw

There is a Maple command to compute an interpolating polynomial:

''f''~((V:= [6,7,-1])) =~ map[3](CurveFitting:-PolynomialInterpolation, [$1..5], 1/~[$1..5], V);

This uses Maple 2018 syntax.

Regarding your handle, "lovemaple": You should show your love by writing some Maple code.

To avoid having a series of Questions with the same title (which is confusing), I changed you titles.

There is a typo in your Poisson probability formula. It should be P(k) = lambda^k*exp(-lambda)/k!.

Both parts of your problem can be done efficiently with a single 4-line procedure that has a combined for-while-loop (thus combining the exercise's requirements for each part), because both parts require summing the above formula for k = 0, 1, 2, ..., N for either a given N (part a) or for an N to be determined by the sum exceeding a certain value (part b). The efficiency is due to not computing lambda^k or k! explicitly; rather, their quotient is built termwise. 

I urge you to deconstruct the 4 main lines (i.e., you can totally ignore the lines that I commented (in the code) as being for advanced users) of the main procedure below, to figure out what each part does (from help pages or asking here), and to reconstruct your own solution from your new-found knowledge. You may also need a brief tutorial or refresher on discrete probability distributions. Wikipedia is an excellent source, and you can also ask here, even if the question has little-to-none Maple content.
  
This code uses Maple 2019 syntax!
 

Download Poisson.mw

In the worksheet below, your 4th proposed solution violates your 5th equation for any m other than -1, 0, or 1. By direct substitution using m=2, p=-1, q=1, it is explicitly shown to violate the 5th equation.

Your 3rd proposed solution is a specialization (to your proposed value of a1) of the 2nd solution returned by solve, which leaves a1 arbitrary.

Your 1st and 2nd proposed solutions, if you rationalize their denominators, can each be seen to be equivalent to a solution returned by solve.
 

Download RatDenom.mw

It's very difficult to measure the time of such trivial operations because the time of the measuring and looping code overwhelms what's being measured.

My first guess is that nops(data) <> 0 is the fastest way. My close second guess would be

not (data::[]).

Here are some hints:

1. There are four ways that may occur to you:

1. Generate the fibonaccis and test if they're square.
2. Generate the squares and test if they're fibonacci.
3. Generate both and find the intersection.
4. Find a formula that just gives square fibonaccis.

Of these, 1 and 3 are easy, and 1 is the most computationally efficient. 2 is more difficult because it's much easier to check if a number is square than to check if it's fibonacci, and also because there are far more squares than fibonaccis. 4, if it were possible at all, would require some sophisticated mathematical work.

2. The predicate issqr tests whether an integer n is a square, i.e., issqr(n) returns true or false depending on whether n is square.

3. A while-loop can be used to stop generation when you exceed 10000.

4. In Maple 2019, the whole thing can be done in a single line of code.