Here is a simple example of the equation with the explicit solution, in which the same error:

dsolve({diff(y(x),x)*(x-1)^2+1=0, y(0)=0});

sol:=dsolve({diff(y(x),x)*(x-1)^2+1=0, y(0)=0}, numeric);

sol(2);

 It is obvious that the curve going from point (0,0) can not be continuously extended further right of  x = 1

plottools:-getdata  command seems appeared in Maple 15. For older versions instead of  plottools:-getdata(A)[3];  you can write

op([1, 1], A);

First, we find 3 specific points on the surface, then - the equation of the plane through these points, and then we prove that each point of the surface lies on this plane:

restart;

r:=[(u-v)^2, u^2-3*v^2, (1/2)*v*(u-2*v)]:

P1, P2, P3:=seq(eval(r,p), p=[[u=0,v=0],[u=1,v=0],[u=0,v=1]]);

LinearAlgebra[Determinant](Matrix([[x,y,z], P2, P3]))=0;

expand(eval(%,[x,y,z]=~r));

Your system is inconsistent for any  rho<>0 . First we eliminate  rho  from equations 2 and 3 and then solve the system of two equations:

restart;

sys:=[8*g(t)^3*diff(g(t),t$2)+4*(g(t)*diff(g(t),t))^2+1=0,rho=-1/g(t)-2*(diff(g(t),t)+t*diff(g(t),t$2))-t/(2*g(t)^3),rho=(-t/g(t))*(diff(g(t),t))^2+t/(4*g(t)^3)];

R:=eliminate(sys[2..3], rho);

Sol:=dsolve([sys[1],op(R[2])]);

eval(sys,Sol[1,1]);

eval(sys,Sol[1,2]);

 For  rho=0  solutions are  Sol .

In the procedure  linalg  package was used but it is not called. Add the line  uses linalg; after local-line of your procedure:

Cocycle:=proc(L,n)

local i,j,k,h,v,u,w,C,eqns,e,f,g;  

uses linalg;

... 

nm

At first I misunderstood your demands. The new procedure  P1  solves the problem for all selected variables in any powers (including the symbolic powers). Note that  simplify  with sideral relation does not work for symbolic powers.

restart;

P1:=proc(expr, S::list)

local Rule;

Rule:=product(S[i]^p[i]::anything, i=1..nops(S))=0;

applyrule(Rule, normal(expand(expr)));

end proc: 

Examples of use:

f0 := a^4+4*a^3*b+6*a^2*b^2+4*a*b^3+b^4:

f1 := 3*(a*b-2*c):

f2 := (a*b-2*c)/(-a*b+c):

f3:=a^b*c:

f4:=a+b*c:

f5:=c+a^m*b^n:

P1(f0,[a,b]); P1(f1,[a,b]); P1(f2,[a,b]); P1(f3,[a,b]); P1(f4,[a,b]); P1(f5,[a,b]);

restart;

P := expr->evalindets(normal(expand(expr)), And(`*`, satisfies(t->nops(indets(t))>1)), t->0): 

Examples of use:

f0 := a^4+4*a^3*b+6*a^2*b^2+4*a*b^3+b^4:

f1 := 3*(a*b-2*c):

f2 := (a*b-2*c)/(-a*b+c):

P(f0);  P(f1);  P(f2);

This works for any rational expressions in several variables, if the denominator is not equal to 0 after substitutions.

In your code you are using a point instead of a multiplication sign many times. For example compare the outputs:

simplify(sin(t).sin(t)+cos(t).cos(t));

simplify(sin(t)*sin(t)+cos(t)*cos(t));

In Maple a point is used for the multiplication of matrices and vectors, e.g.

V:=<a, b, c>:

V.V  assuming real;

If you replace  arg  by  args   the procedure will work but improperly. Error reason: in the body of the procedure you are referring to the argument of the procedure, which is unchanged. See

test([y, y*z-x, -15*x*y-x*z-x], [x,y,z]);

test([y, y*z-x, -15*x*y-x*z-x], [1,1,1]);

                                             [y, y*z-x, -15*x*y-x*z-x]

                                             [y, y*z-x, -15*x*y-x*z-x]

The following version works properly:

test1:=proc()

local f;

f:=args[1]; # a list

[f[2],f[2]*f[3]-f[1],-15*f[1]*f[2]-f[1]*f[3]-f[1]];

end:

or (closer to your version):

restart;

test2:=proc()

local f, vars;

f:=args[1]; # a list

vars:=args[2]: # a list

eval(f, [x,y,z]=~vars);

end:

Examples:

test1([x,y,z]);

test1([1,1,1]);

                               [y, y*z-x, -15*x*y-x*z-x]

                                          [1, 0, -17]

Edited.

May be you want to substitute  (-2+d1), (-1+d1)  terms by new symbols:

subs({d1-1=a, d1-2=b}, y);

nb:=(k-sigma+1)*lambda*L*(gamma*upsilon-delta__1122^k*tau)*upsilon*tau*v/(f__F*sigma*(-tau^2+upsilon^2)*k);

hh := (L*lambda*(k-sigma+1)*upsilon*tau*v)/(f__F*sigma*(-tau^2+upsilon^2)*k)=rho;

simplify(subs(L=solve(hh, L), nb));

S:=a^4*b*c*x^2+a^4*b*c*y+a^4*b*c;

subs(t=a^4*b*c, collect(algsubs(a^4*b*c=t, S), t));

Use LinearAlgebra[DotProduct]  with the option  conjugate=false

Addition:  if you still use linalg[dotprod]  then use  'orthogonal'  option.

Example:

linalg[dotprod](vector([a,b,c]), vector([x,y,z]), 'orthogonal');

                                        a*x+b*y+c*z

Edited.

If a family of curves corresponds to the differential equation  F(x, y, y') = 0 , the orthogonal family will correspond to equation   F(x, y, -1/y') = 0 .

Example:

Equations  y'=2*x  and  -1/y'=2*x .  Their solutions  are  y=x^2+C  and  y=-1/2*ln(abs(x))+C .

Their plots:

A:=plot([seq(x^2+C, C=0..10)], x=-2..2, color=red):

B:=plot([seq(-1/2*ln(abs(x))+C, C=0..10)], x=-2..2, color=blue):

plots[display](A,B, scaling=constrained, view=[-2..2,-1..7]);

A more naive method:

e:= g^((2*(-sigma+k+1))/(-1+sigma))-tau^2:

ex:=(-sigma+k+1)/(-1+sigma):

subs(a=ex, factor(expand(eval(e, ex=a))));

Another way (without  ex):

e:= g^((2*(-sigma+k+1))/(-1+sigma))-tau^2:

factor(applyrule(g^(2*a::anything)=(g^a)^2, e));

Edited.