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

For numerical solutions see  ?plots[odeplot]  and  ?DEtools[DEplot] . If you got the solution explicitly, ie. {x(t)=..., y(t)=..., ...} then  assign(%)  and then  plot([x(t), y(t)],...)  .

I do not think that  PlotPositionVector can solve your problem. A special procedure is required.
The procedure  CurvVector  returns  the list of coordinates the normalized curvature vector for a curve defined by the list of coordinates  L :

restart;

CurvVector:=proc(L)

local e1, e1n, e2, e2n, t;

uses Student[VectorCalculus], LinearAlgebra;

t:=indets(L)[1];

e1:=diff(L,t):

e1n:=unapply(simplify(expand(e1/Norm(convert(e1, Vector),2))), t) assuming t::realcons;

e2:=expand(diff(L,t,t)-(convert(diff(L,t,t),Vector).convert(e1n(t), Vector))*e1n(t));

e2n:=unapply(simplify(expand(e2/Norm(convert(e2, Vector),2))), t) assuming t::realcons;

end proc;

Your example:

with(Student[VectorCalculus]):

L:=[p*cos(p^2), p*sin(p^2)]:

f:=CurvVector(L);

V:= RootedVector(root=eval(L, p=1.5), -f(1.5));

A:=PlotVector(V, color=green):

B:=plot([p*cos(p^2), p*sin(p^2), p=1..2], thickness=2):

plots[display](A, B, scaling=constrained);

restart;

f:=x->tan(x)-x-1:  g:=D(f):

x[0]:=1.5: x[1]:=x[0]-f(x[0])/g(x[0]): epsilon:=10^(-3):

for n while evalf(abs(x[n]-x[n-1]))>=epsilon do

x[n+1]:=x[n]-f(x[n])/g(x[n])

od:

n;  # number of iterations

x[n];  # the root

                                                    7

                                          1.132267726

Verification:

fsolve(tan(x)-x-1);

                                           1.132267725

Actual accuracy is 10^(-9)

restart;

P := x->x^7+14*x^4+35*x^3+14*x^2+7*x+88;

x0 := (1+sqrt(2))^(1/7)-(-1+sqrt(2))^(1/7)-(3+2*sqrt(2))^(1/7)-(3-2*sqrt(2))^(1/7);

is(P(x0) = 0);

factor(P(x), (1+sqrt(2))^(1/7)):

x := solve(op(3, %));

is(x = x0);

Your system can be solved exactly by  dsolve  command with  method = laplace  option. The solutions obtained are expressed through the roots of the fifth degree polynom. Then they calculated approximately and simplified:

res := dsolve(`union`({EQ1, EQ2, EQ3, EQ4, EQ5}, {q1(0) = 0.1e-2, q2(0) = 0, q3(0) = 0, q4(0) = 0, q5(0) = 0, (D(q1))(0) = 0, (D(q2))(0) = 0, (D(q3))(0) = 0}), {seq((q || i)(t), i = 1 .. 5)}, method = laplace);

sol := evalf(res);

evalc(simplify(fnormal(sol, 7), zero));

assign(%);

You have a linear system of differential equations in matrix form.

Example of solving:

Y(t):=Vector([u(t),v(t)]):

A:=Matrix([[1,2],[3,4]]):

dsolve({diff(Y(t)[1],t)=(A.Y(t))[1],diff(Y(t)[2],t)=(A.Y(t))[2]} );

You specify an undirected graph, so the distance matrix should be symmetric. But in your graph, for example,  (2, 22) = 547  and  (22, 2) = 500