I want to generate a system of differential equations for DSolve. The unknown function p[t] is a 2D-function with parameter t and coordinates {x[t], y[t]}. Think of the movement of a point in a plane.

Thereby I encountered a general problem: As long as p is not represented by a list, I cannot tell WL that p is a vector, not a scalar.

I think there should be a symbol attribute, telling that the symbol is meant to be a nD-vector (or nxm matrice, etc.).

See the following piece of code, where p0={p0x, p0y} is a fixed vector, and p shall be a general 2D-Vector. I want to calculate and use d(p) = Norm(p0 - p) within a larger program:

Block[{p0, p0x = 11, p0y = 12, p, x, y, d},
  p0 = {p0x, p0y};
  Print["p0 = ", p0, ", Attr(p0) = ", Attributes[p0], "\tp = ", p, 
   ", Attr(p) = ", Attributes[p]];
  d = Norm[p0 - p];
  Print["1. d(p) = Norm({p0x,p0y} - p) = ", d];
  p = {x, y}; Print["2. d(p=", p, ") = ", d,
   {x, y} = {21, 22}; "\td(p=", p, ") = ", d];
  Clear[x, y];
  d = Norm[p0 - p];
  Print["3. d(p=", p, ") = Norm({p0x,p0y} - {x,y}) = ", d,
   {x, y} = {21, 22}; "\td(p=", p, ") = ", d];
  ];

p0 = {11,12}, Attr(p0) = {} p = p, Attr(p) = {}

1. d(p) = Norm({p0x,p0y} - p) = Sqrt[Abs[11-p]^2+Abs[12-p]^2]

2. d(p={x,y}) = {Sqrt[Abs[11-x]^2+Abs[12-x]^2],Sqrt[Abs[11-y]^2+Abs[12-y]^2]}   d(p={21,22}) = {Sqrt[181],Sqrt[221]}

3. d(p={x,y}) = Norm({p0x,p0y} - {x,y}) = Sqrt[Abs[11-x]^2+Abs[12-y]^2] d(p={21,22}) = 10 Sqrt[2]

In the first step, I calculate Norm[p0 - p], which gives the correct result for a scalar p, but a wrong result for a vector p. If in the second step I replace p by {x,y} within the result of the first step, I get the wrong result of course.
In the third step, I recalculate Norm[p0 - p], which is correct now.

Of course I can redefine d(p) by a function with set-delayed: d[p_] := Norm[p0 - p] Then everything works fine, if I use d[p] instead of Norm[p0 -p] and replacement.

Block[{p0, p0x = 11, p0y = 12, p, x, y, d},
  p0 = {p0x, p0y};
  Print["p0 = ", p0, ", Attr(p0) = ", Attributes[p0], "\tp = ", p, 
   ", Attr(p) = ", Attributes[p]];
  d[p_] := Norm[p0 - p];
  Print["1. d(p) = Norm({p0x,p0y} - p) = ", d[p]];
  p = {x, y}; Print["2. d(p=", p, ") = Norm({p0x,p0y} - p] = ", d[p],
   {x, y} = {21, 22}; "\td(p=", p, ") = ", d[p]];
  Clear[x, y];
  Print["3. d(p=", p, ") = Norm({p0x,p0y} - {x,y}) = ", d[p],
   {x, y} = {21, 22}; "\td(p=", p, ") = ", d[p]];
  ];

p0 = {11,12}, Attr(p0) = {} p = p, Attr(p) = {}

1. d(p) = Norm({p0x,p0y} - p) = Sqrt[Abs[11-p]^2+Abs[12-p]^2]

2. d(p={x,y}) = Norm({p0x,p0y} - p] = Sqrt[Abs[11-x]^2+Abs[12-y]^2] d(p={21,22}) = 10 Sqrt[2]

3. d(p={x,y}) = Norm({p0x,p0y} - {x,y}) = Sqrt[Abs[11-x]^2+Abs[12-y]^2] d(p={21,22}) = 10 Sqrt[2]

The problem is, that this doesn't help in my original problem of a differential equation for a vector-valued function p[t].

See the following code for a trivial constant-velocity movement in a plane, starting at some point p0:

Block[{p0 = {1, 2}, t, p, v = {2, 1}, eqs, sol},
  eqs = {p'[t] == v, p[0] == p0};
  Print["eqs = ", Column[eqs]];
  sol = DSolve[eqs, p[t], {t, 0, 4}];
  Print["sol = ", sol];
  ];

eqs = (p^\[Prime])[t]=={2,1}
p[0]=={1,2}

DSolve::nolist: List encountered within {(p^\[Prime])[t]=={2,1}}. There should be no lists on either side of the equations.

sol = DSolve[{(p^\[Prime])[t]=={2,1},p[0]=={1,2}},p[t],{t,0,4}]

Even though by the boundary condition p[0]==p0 it is clear that p is a vector, WL doesn't recognize p[t] and therefore p'[t] to be vector-valued functions of a scalar t.

Since I cannot tell WL that p[t] is vector-valued, I am forced to switch to coordinate-based differential equations:

Block[{p0 = {1, 2}, t, p, x, y, v = {2, 1}, eqs, sol},
  p[t_] = {x[t], y[t]};
  eqs = {p'[t] == v, p[0] == p0};
  Print["eqs = ", Column[eqs]];
  sol = DSolve[eqs, p[t], {t, 0, 4}];
  Print["sol = ", sol];
  ];

eqs = {(x^\[Prime])[t],(y^\[Prime])[t]}=={2,1}
{x[0],y[0]}=={1,2}

sol = {{x[t]->1+2 t,y[t]->2+t}}

This seems to be no big issue but imagine p'[t] as a complicated expression in p[t] and t. You have to switch to NDSolve then and you will get a list of scalar InterpolatingFunction(s) as the solution for {x[t], y[t]}. It would be much better to get one vector-valued InterpolationFunction as the solution for p[t].