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

You can easily find the general solution of this equation:

restart;

g :=0.88641:

e := 2.53128:

eq := tan(g)-e*sin(f)/(1+e*cos(f)):

solve(eq, AllSolutions);

                    -2.566269006+6.283185307*_Z1,   1.197496352+6.283185307*_Z1

about(_Z1);

                                       Originally _Z1, renamed _Z1~:
                                          is assumed to be: integer

restart;

a := 0:  b := 2:

n[x] := 20:  n[t] := 50:  dt := 0.01:  c := 1:  dx := (b-a)/n[x]:

for i while i <= n[x] do

  if 5 <= i and i <= 10 then u[i] := 2 else u[i] := 1:

  end if;

 end do;

seq(u[i], i = 1 .. 20);

                                    1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1

P:=proc(N::posint, F::procedure)

uses plottools;

seq(seq(disk([i,j], 0.1, color=`if`(F(i,j), red, blue)), i=1..N), j=1..N);

end:

Example:

plots[display](P(7, (i, j)->`if`(i>j, true, false)), axes=none);

 Improvement.

 Instead

plots[display](P(7, (i, j)->`if`(i>j, true, false)), axes=none); 

it is shorter to write 

plots[display](P(7, (i, j)->is(i>j)), axes=none);