Your functions  b[n,m]  is easy to construct by a double loop. Procedure  HybrFunc  solves this problem. Functions  b[n,m]  can be called by their names (b - global variable) or as array elements. You have not written what values variable  tj  takes , so I just specify it how procedure argument.

restart;

HybrFunc:=proc(N, M, tj)

local T, n, m;

global b;

for n to N do

for m from 0 to M-1 do

T[m]:=t->t^m;

b[n,m]:=unapply(piecewise(t>=(n-1)*tj/N and t<n*tj/N, T[m](N*t-(n-1)*tj), 0), t);

od; od;

Array(1..N, 0..M-1, (n,m)->b[n,m](t));

end proc:

Example:

HybrFunc(3, 4, 10)[3,1];

b[3,2](t);

restart;

sort(x+y+z, [x, y, z], ascending);

                     z+y+x

restart: 

Mf(x):=piecewise(x<=L/2,1/2*x*F,x>1/2*L,1/2*x*F-F*(x-1/2*L)): 

eq[1]:=Mf(xi)*F*L=Mf(x): 

Mf(xi):=simplify(convert(solve(eq[1],Mf(xi)), piecewise,x)); 

Mf(xi):=subs(x=xi*L,Mf(xi)); 

Mf(xi):=simplify(Mf(xi)) assuming F>0, L>0;

plot(x^2, x=-1..2, tickmarks=[0,0]);

f := int((A*sin(omega*t+phi)-B*sin(omega*t))^2, t):

g := int((A*sin(omega*t)*cos(phi)+A*cos(omega*t)*sin(phi)-B*sin(omega*t))^2, t):

delta=simplify(expand(f)-expand(g));  # Difference between f and g 

We see that the difference between  f  and  g  does not depend on  t, so there is no contradiction.

ex := taylor(exp(x), x = 0, 5);

eval(ex, x=1);

subs(x=1, convert(ex, polynom));

A:=Matrix(3,2, symbol=a):

B:=Matrix(2,3, symbol=b):

C:=Matrix([[8, 2, -2], [2, 5, 4], [-2, 4, 5]]):

simplify(B.A, {seq(seq((A.B)[i,j]=C[i,j], j=1..3), i=1..3)});

mahmood180, I tried to upload your file, but the 404 error appeared.

Look at my solution:

restart;

SpecialMatrix := proc (A::list)

local r, a, b, A1, A2, a1, a2; 

r := (1/2)*nops(A)-1/2;

assign(seq(a[i-1] = A[i], i = 1 .. r+1), seq(b[i] = A[r+1+i], i = 1 .. r));

A1[1, 1] := Matrix(1, {(1, 1) = 2*a[0]});

A1[1, 2] := Matrix([[seq(a[i], i = 1 .. r)]]);

A1[1, 3] := Matrix([[seq(b[i], i = 1 .. r)]]); A1[2, 1] := A1[1, 2]^%T;

A1[2, 2] := Matrix(r, {seq(seq((i, j) = a[j-i], j = i+1 .. r), i = 1 .. r-1), seq((i, i) = 2*a[0], i = 1 .. r)}, shape = symmetric);

A1[2, 3] := Matrix(r, {seq(seq((i, j) = b[j-i], j = i+1 .. r), i = 1 .. r-1), seq((i, i) = 0, i = 1 .. r)}, shape = antisymmetric);

A1[3, 1] := A1[1, 3]^%T; A1[3, 2] := -A1[2, 3]; A1[3, 3] := A1[2, 2];

a1 := `<,>`(seq(`<|>`(A1[i, 1], A1[i, 2], A1[i, 3]), i = 1 .. 3)); 

A2[1, 1] := Matrix(1);

A2[1, 2] := A1[1, 2];

A2[1, 3] := A1[1, 3];

A2[2, 1] := Matrix(r, 1);

A2[2, 2] := Matrix(r,{seq(seq((i, j) = a[i+j], j = 1 .. r-i), i = 1 .. r-1)});

A2[2, 3] := Matrix(r,{seq(seq((i, j) = b[i+j], j = 1 .. r-i), i = 1 .. r-1)});

A2[3, 1] := A2[2, 1]; A2[3, 2] := A2[2, 3]; A2[3, 3] := -A2[2, 2];

a2 := `<,>`(seq(`<|>`(A2[i, 1], A2[i, 2], A2[i, 3]), i = 1 .. 3)); 

(1/2)*a1+(1/2)*a2; 

end proc:

Example for r=4:

SpecialMatrix([seq(a[i], i = 0 .. 4), seq(a[i]^`&lowast;`, i = 1 .. 4)]);

SymmetricSum:=proc(S::{list,set}, P::list)

local L;

uses combinat;

if nops(S)<>nops(P) then error "Should be nops(S)=nops(P)" fi;

L:=permute(P);

add(mul(S[i]^l[i], i=1..nops(S)), l=L);

end proc;

CyclicSum:=proc(S::list, P::list)

local n, L0, L;

if nops(S)<>nops(P) then error "Should be nops(S)=nops(P)" fi;

n:=nops(P);

L0:=[op(P),op(P)]; L:=[seq([seq(L0[i+j], j=0..n-1)],i=1..n)];

add(mul(S[i]^l[i], i=1..nops(S)), l=L);

end proc:

Examples:

SymmetricSum([x,y,z], [2,3,5]);

CyclicSum([x,y,z,u], [2,3,4,5]);

Formal parameters:  L - the list of your functions  psi ,  m - a symbol or a number,  M - a positive integer.

MMatrix := proc (L::list, m, M::posint)

Matrix(M, [seq([seq(L[i]((1/2)*(2*j-1)/m), j = 1 .. M)], i = 1 .. M)]);

end proc:

Example:

MMatrix([seq(psi[1, i], i = 0 .. 4)], 10, 5);

Of course instead of arbitrary functions psi, you can write as the entries of the list L your specific functions.

Another way of filling the interior of the cardioid used. It is based on an approximation by polygons. The same method was used in my procedure Picture. See  http://www.mapleprimes.com/posts/145922-Perimeter-Area-And-Visualization-Of-A-Plane-Figure-

restart;

N := 120:

r := 2-2*sin(2*Pi*k/N):

L := [[0, 0], seq([r*cos(2*Pi*k/N), r*sin(2*Pi*k/N)], k = 0 .. N)]:

A := seq(plottools[polygon](L[1 .. n], color = cyan), n = 2 .. N+2):

B := seq(plottools[curve](L[2 .. n], color = black, thickness = 3), n = 2 .. N+2):

l := plottools[line]([0, 0], [2, 0], color = cyan, thickness = 4):

C := seq(plots[display](l, A[i], B[i]), i = 1 .. N+1):

plots[display](C, insequence = true, axiscoordinates = polar, scaling = constrained);

tmp:=[[0, 1, 2], [1, 0, 2], [1, 1, 2], [1, 2, 0]];

select(x->is(`+`(op(x))=3), tmp);

                            tmp := [[0, 1, 2], [1, 0, 2], [1, 1, 2], [1, 2, 0]]

                                       [[0, 1, 2], [1, 0, 2], [1, 2, 0]]

triType:=proc(a,b,c)

local A, B, C;

if is(a+b<=c) or is(a+c<=b) or is(c+b<=a) then error "No triangle" fi;

A:=b^2+c^2-a^2; B:=a^2+c^2-b^2; C:=a^2+b^2-c^2;

if is(A>0) and is(B>0) and is(C>0) then return acute else

if is(A=0) or is(B=0) or is(C=0) then return right else

obtuse fi; fi;

end proc;

Examples:

triType(3, 3, 4);

triType(3, 4, 5);

triType(3, 4, 6);

triType(3, 5, 10);

                                  acute

                                   right

                                 obtuse

Error, (in triType) No triangle

restart;

B:= ln(exp(t)+a+sqrt(exp(t)-a+b)):

e:=solve(A(t)=B, exp(t)):

subs(exp(t)=e[1], diff(B, t));

simplify(x^4+y^2+2*y*x^2, {x^2+y=z});

                                  z^2

or

algsubs(x^2+y=z, x^4+y^2+2*y*x^2);

                                  z^2

or

compoly(x^4+y^2+2*y*x^2, {x, y});   # composition of polynomials

                        x^2, x = x^2+y