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);

plots[implicitplot]([seq(x^2-2*y^2=C, C=-5..5)], x= -5..5, y= -5..5, numpoints=10000, scaling=constrained);

Numpoints  option allows you to significantly improve the quality of construction, and  scaling=constrained  option equalizes the scale of the axes.

Eq1:=diff(f(x), x)=1, g(x)=diff(f(x), x);

cs:=f(0)=0;

p:=dsolve({Eq1, cs}, numeric);

plots[odeplot](p, [x, g(x)], 0..10);

If you open the brackets in the canonical equation of a circle  (x-x0)^2+(y-y0)^2=R^2 , you get an equation of the form  x^2+y^2+a*x+b*y+c=0 . 

Your problem can be easily solved symbolically, ie exactly. Here is an example:

restart;

F := x^2+y^2+a*x+b*y+c:

XY := [[0, 0], [2, 3/2], [3, 1], [2, -1]]:

minimize(add((eval(F, {x = XY[i, 1], y = XY[i, 2]}))^2, i = 1 .. nops(XY)), location);

assign(op(%[2])[1]):

A := plot(XY, style = point, symbolsize = 15, symbol = solidcircle, color = red):

B := plots[implicitplot](F, x = -1 .. 4, y = -2 .. 2, color = green):

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

Parameters  x0, y0, R  can be easily expressed in terms of  a, b, c . This you can make yourself.

plots[implicitplot](0.3939708949*x^2 - 0.005975799853 + 0.6345432059*y^2=0,

x=-1..1, y=-1..1, thickness=2, numpoints=50000, scaling=constrained);

Addition.  We can do the job without  numpoints option if to shorten the axis, because the ellipse is very small:

plots[implicitplot](0.3939708949*x^2 - 0.005975799853 + 0.6345432059*y^2=0,

x=-0.2..0.2, y=-0.2..0.2, thickness=2, scaling=constrained);