If I understand your question, you want the maximum element in the matrix put in the needed position  and use as a pivoting element. This can be done by interchanging of rows and/or columns. The procedure  MaxPivot sets maximum element of matrix  A  at the position, defined by the list  L:

MaxPivot:=proc(A::Matrix, L::list)

local m, n, M, i, j;

uses LinearAlgebra;

m:=RowDimension(A); n:=ColumnDimension(A);

M:=[[0,0], 0];

for i to m do

for j to n do

if abs(A[i,j])>M[2] then M:=[[i,j], abs(A[i,j])] fi;

od; od;

RowOperation(A, [L[1], M[1,1]]);

ColumnOperation(%, [L[2], M[1,2]]);

end:

Example:

A:=<<1 | 2 | 4> , <5 | 6 | 9> , <-1 | 3 | 4>>:
MaxPivot(A, [1,1]);

The basic idea: if  a<b  that of the inequality  f(x)<g(x)  imply that  int(f(x), x=a..b)<int(g(x), x=a..b) . Sufficient to prove the  original limit for x> 0 .

For  0<x<Pi  consider the obvious double inequality  0<sin(x)<x  and integrate  it several times from 0 to x , expressing everytime  sin(x)  or  cos(x) . This can easily be automated in Maple:

L:=[0<sin(x), sin(x)<x]:

for k to 4 do

L:=[seq(op(solve(map(int, L[i], x=0..x), `if`(k::odd, cos(x), sin(x)))), i=1..2)];

od;  

From the last inequality we express  A = (x-sin (x)) / x^3 :

seq(solve(subs(sin(x)=x-A*x^3,L[i]), {A}) , i=1..2) assuming x^3>0;

So we have a double eniquality  -1/120*x^2+1/6 < A <1/6  , from which it obviously follows that

limit(A, x=0)=1/6

int(sqrt(diff(sqrt(20)*cos(t),t)^2+diff(sqrt(13)*sin(t),t)^2), t=0..2*Pi);
evalf(%);

                                     8*5^(1/2)*EllipticE(1/10*35^(1/2))
                                                      25.44987267

Did you remember to put a semicolon or a colon at the end of the first line?

NumberOfDigits:=proc(N)

local X, n;

X[0]:=0; X[1]:=1; X[2]:=2;

for n from 3 to N do

X[n]:=n*(X[n-3]+X[n-2]+X[n-1]);

od;

convert(X[N], string);

length(%);

end proc;

Your example:

NumberOfDigits(2013);

                  5782

Position:=proc(x::anything, L::list)

local p;

if member(x, L, 'p') then p else 0 fi;

end proc:

Examples:

Position(2, [3,8,2,7]);

Position(5, [3,8,2,7]);

                    3

                    0

minimize(K(N), N=1..5, location);

maximize(K(N), N=0.5..1, location);

                3.019929556, {[{N = 1.279315996}, 3.019929556]}
                3.663102178, {[{N = 0 .6129946770}, 3.663102178]}

F:=0.85:  B:=0.5:

K:=N->N*(1+F*N/(N^2+B^2-F*N));

plot(K, 0..10, thickness=2);

minimize(K(N), N=0..10);

maximize(K(N), N=0..10);

minimize(K(N), N=1..5);

maximize(K(N), N=0.5..1);

The procedure  Fourier  solves your problem. Formal arguments:  f  is original function,  if you need the whole series then  n  is  infinity ,  if the partial sum then  n  is nonnegative integer.

Fourier:=proc(f, n::{nonnegint, infinity})

local a0, a, b;

a0:=(1/(2*Pi))*(int(f(x),x=-Pi..Pi));

a:=k->(1/Pi)*(int(f(x)*cos(k*x),x=-Pi..Pi, AllSolutions) assuming k::posint);

b:=k->(1/Pi)*(int(f(x)*sin(k*x),x=-Pi..Pi, AllSolutions) assuming k::posint);

if n=infinity then return a0+Sum(simplify(a(k)*cos(k*x)+b(k)*sin(k*x), assume=k::posint) , k=1..n) else a0+add(a(k)*cos(k*x)+b(k)*sin(k*x), k=1..n)  fi;

end proc;

Examples:

f := x->min(abs(x), (1/2)*Pi);

Fourier(f, infinity);

Fourier(f, 5);

plot([f(x), seq(Fourier(f, k), k = [1, 2, 3, 5, 6, 30])], x = -Pi .. Pi, thickness = [3, `$`(1, 6)], color =

 [blue, red, brown, khaki, yellow, green, violet], scaling = constrained);

 PS.  For the violet  curve (k = 30)  the plot is indistinguishable from the original one.

Edited: The code is corrected.

A:=1/2*pochhammer(1-n,n):

eval(A, n=0), eval(A, n=k) assuming k::posint;

                                  1/2,  0

plots[implicitplot](y*(1-x-y)>=0, x=0..1, y=0..1,  coloring=[green, white], filledregions, numpoints=10000, view=[-0.5..1.5, -0.5..1.5]);

sort([[2,1],[3,2],[5,3],[1,4],[4,5]], (a,b)->is(a[1]<=b[1]));

                     [[1, 4], [2, 1], [3, 2], [4, 5], [5, 3]]

Replace the penultimate line of your code by

sort([seq(B[i], i=1..3)], proc (x, y) options operator, arrow; is(y[2] < x[2]) end proc);

L:=[[0,2,3,0,0], [2,3,0,0,0], [1,1,2,2,0], [2,2,1]]:

convert(map(x->1/mul(x[i]!, i=1..nops(x)), L), `+`);

                                        2/3

B1:=[2., 36.718220544331125]:

B2:=[2., 35.08378362904457]:

B3:=[2., 51.78712780854305]:

sort([B1, B2, B3], (x,y)->is(x[2]>y[2]));

datamax:=%[1];