Good question.  It is a habit from using the older matrix structure.  There eval was needed because matrices (small-m), as well as tables, arrays, procedures, and modules (records) use last-name-evaluation, which means that if M were a (small-m) matrix, then returning M would actually return the local name M rather than the matrix M.  Because rtables (of which Matrix is a specific type) do not use last-name-evaluation this isn't necessary. 

See the help page ?last_name_eval for details.

FYI, your procedure had a couple of flaws.  First, it never assigned a Matrix to M.  Second, it never returned the computed Matrix.  Here is a corrected version.  Of course, the methods that have been suggested are more convenient.

mat:= proc(a,b,c,N)
local i,j,M:
    M := Matrix(N,N);
    for i from 1 to N do
        for j from 1 to N do
            if j>(i+1) then M[i,j]:= a:
            elif i>(j+1) then M[i,j]:= c:
            else M[i,j]:= b:
            end if:
        end do:
    end do:
    eval(M);
end proc:
mat(a,b,c,4);

> Matrix(4,4,(i,j)->`if`(j>i+1,a,`if`(i>j+1,c,b)));   
                                  [b    b    a    a]
                                  [                ]
                                  [b    b    b    a]
                                  [                ]
                                  [c    b    b    b]
                                  [                ]
                                  [c    c    b    b]

eq1 := x^3 + y^3 = 3*x^2*y:
eq2 := subs(y=y(x),eq1);
                                 3       3      2
                         eq2 := x  + y(x)  = 3 x  y(x)

isolate(diff(eq2,x),diff(y(x),x));
                                        2
                          d         -3 x  + 6 x y(x)
                          -- y(x) = ----------------
                          dx               2      2
                                     3 y(x)  - 3 x

subs(y(x)=y,rhs(%));
                                     2
                                 -3 x  + 6 x y
                                 -------------
                                     2      2
                                  3 y  - 3 x

alias(y=y(x)):
eq := x^3 + y^3 = 3*x^2*y:
isolate(diff(eq,x),diff(y,x));

                                        2
                             d      -3 x  + 6 x y
                             -- y = -------------
                             dx         2      2
                                     3 y  - 3 x

First off, your definition of coprime is nonstandard (putting it politely).

I'd try

iscoprime := (k,j) -> evalb(igcd(k,j)=1):
L := [seq(1..65)]:
select(iscoprime, L, 66):

If you suspect that the answer is a rational polynomial, you can try to solve for the coefficients

f := 1/(k^2+3*k+1):

nn := 1:
nd := 3:

g := add(a||i*k^i,i=0..nn)/add(b||i*k^i,i=0..nd);
                                     a0 + a1 k
                        g := -------------------------
                                             2       3
                             b0 + b1 k + b2 k  + b3 k


eqs := {seq(numer(eval(g-f,k=kk)),kk=1..nn+nd+2)}:
sol := solve(eqs);
sol := {a0 = -3 b3 + b2, a1 = b3, b0 = -3 b3 + b2, b1 = -8 b3 + 3 b2, b2 = b2,

    b3 = b3}

normal(subs(sol,g));
                                      1
                                 ------------
                                  2
                                 k  + 3 k + 1

Here's a useful heuristic.  Given a problem that you cannot solve, create an easier version of it that you can, than attempt to extend the solution.  That is the reason you were asked to solve for p(2) and p(3) before p(n).  However, you can make it even simpler.  The general problem is, given N days in a year ... find p(n), where n is the number of  people ....   Suppose N=1.  What is p(2)?   p(3)?  p(n)?  Now suppose N=2.

That the ring and the sphere have different Gaussian curvatures directly proves that there will be a "bucklement", however, if you replace the sphere with a cone you will also get a "bucklement" even though the cone has zero Gaussian curvature.  You can, however, wrap the ring around the cone (rather than encircling it and then trying to squash it flat).

M := <<1|2>,<3|4>>:
LinearAlgebra:-MatrixPower(M,4);
                                     [199    290]
                                     [          ]
                                     [435    634]

String matching might be the best way (I used 10,000 digits here):

module()
local pi := cat("3",convert(evalf[10\000](Pi),'string')[3..-1]);
    TypeTools:-AddType(inpi
                       , proc(n)
                           n :: 'nonnegint'
                           and (StringTools:-Search(convert(n,'string'),pi) <> 0)
                       end proc
                      );
end module:

type(14122,inpi);
                                true
type(141222,inpi);
                                false

Note that there is a typo in ID4dof.mws:

Test := MatrixInverse(G11-P1):

should be

Test := MatrixInverse(G11-Pt):

also the commented-out line that follows has a typo, "Multply". 

That didn't solve your problem...

Use the unapply procedure to convert an expression into an operator:

with(stats):
does_work:= proc(data)
local d0,d1,x,y;
    unapply(rhs(fit['leastsquare'[[x,y], y=d0+d1*x]](data)),x);
end proc:
f := does_work([[0,1,2],[3,4,5]]);                        
                       f := x -> 3 + x

f(a);
                       3+a

Note that the stats package has been supreceded by the Statistics package.  Using that package you could do

does_work:= proc(data)
local x;
    unapply(Statistics:-LinearFit([1,x],data[],x),x);
end proc:

f := does_work([[0,1,2],[3,4,5]]);
             f := x -> 3.00000000000000044 + 0.99999999999999988 x


f(a);
                  3.00000000000000044 + 0.99999999999999988 a

does_work:= proc(data)
local f,x;
    f :=Statistics:-LinearFit([1,x],data[],x);
    unapply(convert(f,'rational'),x);
end proc:

f := does_work([[0,1,2],[3,4,5]]);
                                f := x -> 3 + x


f(a);
                                     3 + a

Assign L the list of values (1's and 0's). Then do

remove(Testzero, map(nops, [ListTools:-Split(`=`, L, 1)]));

You can do

interface(prettyprint=0):

The downside is that the output is then left aligned, and harder to read. If you only want to convert a particular expression, you can use the procedure lprint.

Weird...my browser won't permit copying the Maple code that you entered.  I never had that problem before.

I don't understand why you are specifying a range from 0 to infinity.  The solution is near 0.