> these complex parts could been multiplied together and give a REAL noise

ah yes that is potentially worrisome, I think. You know I'm totally out of my depth here -- there is a Package called RealDomain, might it help? The help says:

By default, Maple performs computations under the assumption that the underlying number system is the complex field. The RealDomain package provides an environment in which computations are performed under the assumption that the basic underlying number system is the field of real numbers.

and here are two ways to use it:

with( RealDomain );

use RealDomain in simplify( sqrt( x^2 ) ) end use;

> these complex parts could been multiplied together and give a REAL noise

ah yes that is potentially worrisome, I think. You know I'm totally out of my depth here -- there is a Package called RealDomain, might it help? The help says:

By default, Maple performs computations under the assumption that the underlying number system is the complex field. The RealDomain package provides an environment in which computations are performed under the assumption that the basic underlying number system is the field of real numbers.

and here are two ways to use it:

with( RealDomain );

use RealDomain in simplify( sqrt( x^2 ) ) end use;

learned a lot from this discussion, thanks.

fyi, M13 on my machine (Vista) spit out:

[-.47477721457]
[-.996172636207]
[-.5000137071086]
[-.990140076772]

Not sure what you mean.

plus this is a Maple discussion group not LaTeX ...

Not sure what you mean.

plus this is a Maple discussion group not LaTeX ...

Plot associated with:

newsys := [
diff(x(t),t) = (x(t)^(1/2)-1/20*x(t)-c(t))^2, 
diff(c(t),t) = 1/2*c(t)*(1/2*1/(x(t)^(1/2))-1/20-q(t))*(x(t)^(1/2)-1/20*x(t)-c(t)), 
diff(q(t),t) = (-(q(t)-1/50)*(q(t)-3/100)
+1/2*(q(t)-1/40)*c(t)*(1/2*1/(x(t)^(1/2))-1/20-q(t))/(x(t)^(1/2)-1/20*x(t)-c(t)))*(x(t)^(1/2)-1/20*x(t)-c(t))
];

newvars := [
x(t) = (u(t)+400/9)^2, 
c(t) = ((u(t)+400/9)^2)^(1/2)-1/20*(u(t)+400/9)^2-v(t), 
q(t) = 1/2*1/(((u(t)+400/9)^2)^(1/2))-1/20-w(t)/(1/2*((u(t)+400/9)^2)^(1/2)-1/40*(u(t)+400/9)^2-1/2*v(t))
];

Plot associated with:

newsys := [
diff(x(t),t) = (x(t)^(1/2)-1/20*x(t)-c(t))^3, 
diff(c(t),t) = 1/2*c(t)*(1/2*1/(x(t)^(1/2))-1/20-q(t))*(x(t)^(1/2)-1/20*x(t)-c(t))^2, 
diff(q(t),t) = (-(q(t)-1/50)*(q(t)-3/100)
+1/2*(q(t)-1/40)*c(t)*(1/2*1/(x(t)^(1/2))-1/20-q(t))/(x(t)^(1/2)-1/20*x(t)-c(t)))*(x(t)^(1/2)-1/20*x(t)-c(t))^2
];

newvars := [
x(t) = (u(t)+400/9)^2, 
c(t) = ((u(t)+400/9)^2)^(1/2)-1/20*(u(t)+400/9)^2-v(t), 
q(t) = 1/2*1/(((u(t)+400/9)^2)^(1/2))-1/20-w(t)/(1/2*((u(t)+400/9)^2)^(1/2)-1/40*(u(t)+400/9)^2-1/2*v(t))
];

Plot associated with:

newsys := [
diff(x(t),t) = (x(t)^(1/2)-1/20*x(t)-c(t))^2, 
diff(c(t),t) = 1/2*c(t)*(1/2*1/(x(t)^(1/2))-1/20-q(t))*(x(t)^(1/2)-1/20*x(t)-c(t)), 
diff(q(t),t) = (-(q(t)-1/50)*(q(t)-3/100)
+1/2*(q(t)-1/40)*c(t)*(1/2*1/(x(t)^(1/2))-1/20-q(t))/(x(t)^(1/2)-1/20*x(t)-c(t)))*(x(t)^(1/2)-1/20*x(t)-c(t))
];

newvars := [
x(t) = (u(t)+400/9)^2, 
c(t) = ((u(t)+400/9)^2)^(1/2)-1/20*(u(t)+400/9)^2-v(t), 
q(t) = 1/2*1/(((u(t)+400/9)^2)^(1/2))-1/20-w(t)/(1/2*((u(t)+400/9)^2)^(1/2)-1/40*(u(t)+400/9)^2-1/2*v(t))
];

Plot associated with:

newsys := [
diff(x(t),t) = (x(t)^(1/2)-1/20*x(t)-c(t))^3, 
diff(c(t),t) = 1/2*c(t)*(1/2*1/(x(t)^(1/2))-1/20-q(t))*(x(t)^(1/2)-1/20*x(t)-c(t))^2, 
diff(q(t),t) = (-(q(t)-1/50)*(q(t)-3/100)
+1/2*(q(t)-1/40)*c(t)*(1/2*1/(x(t)^(1/2))-1/20-q(t))/(x(t)^(1/2)-1/20*x(t)-c(t)))*(x(t)^(1/2)-1/20*x(t)-c(t))^2
];

newvars := [
x(t) = (u(t)+400/9)^2, 
c(t) = ((u(t)+400/9)^2)^(1/2)-1/20*(u(t)+400/9)^2-v(t), 
q(t) = 1/2*1/(((u(t)+400/9)^2)^(1/2))-1/20-w(t)/(1/2*((u(t)+400/9)^2)^(1/2)-1/40*(u(t)+400/9)^2-1/2*v(t))
];

As I noted, a typo crept up in the system analyzed by Robert. Two days ago I reported my inability to diagonalize the system once the typo had been corrected. I have now found a transformation that diagonalizes the system. I report on it below. I'm still unsure of the interpretation of it all though...

Using a transformation similar to the one Robert applied to the typo-contaminated system (TRANSFORMATION 1 below) yields the non-diagonal matrix:

But the following transformation works:

newsys := [
diff(x(t),t) = (x(t)^(1/2)-1/20*x(t)-c(t))^2, 
diff(c(t),t) = 1/2*c(t)*(1/2*1/(x(t)^(1/2))-1/20-q(t))*(x(t)^(1/2)-1/20*x(t)-c(t)), 
diff(q(t),t) = (-(q(t)-1/50)*(q(t)-3/100)
+1/2*(q(t)-1/40)*c(t)*(1/2*1/(x(t)^(1/2))-1/20-q(t))/(x(t)^(1/2)-1/20*x(t)-c(t)))*(x(t)^(1/2)-1/20*x(t)-c(t))
];

newvars := [
x(t) = (u(t)+400/9)^2, 
c(t) = ((u(t)+400/9)^2)^(1/2)-1/20*(u(t)+400/9)^2-v(t), 
q(t) = 1/2*1/(((u(t)+400/9)^2)^(1/2))-1/20-w(t)/(1/2*((u(t)+400/9)^2)^(1/2)-1/40*(u(t)+400/9)^2-1/2*v(t))
];

Its Jordan form is diagonal:

MATRIX([[0, 0, 0], [0, -187/108, 0], [0, 0, 0]]);

The phase diagram of the 2-D approximation of the transformed system is:

[SERVER UNABLE TO UPLOAD FILES, FIGURE WILL BE UPLOADED LATER BELOW]

My interpretation of the phase diagram is that the stationary state (u=v=w=0) is unstable. Do you agree?

But now consider a different transformation, which allows the derivative of x(t) to be both positive and negative,

newsys := [
diff(x(t),t) = (x(t)^(1/2)-1/20*x(t)-c(t))^3, 
diff(c(t),t) = 1/2*c(t)*(1/2*1/(x(t)^(1/2))-1/20-q(t))*(x(t)^(1/2)-1/20*x(t)-c(t))^2, 
diff(q(t),t) = (-(q(t)-1/50)*(q(t)-3/100)
+1/2*(q(t)-1/40)*c(t)*(1/2*1/(x(t)^(1/2))-1/20-q(t))/(x(t)^(1/2)-1/20*x(t)-c(t)))*(x(t)^(1/2)-1/20*x(t)-c(t))^2
];

newvars := [
x(t) = (u(t)+400/9)^2, 
c(t) = ((u(t)+400/9)^2)^(1/2)-1/20*(u(t)+400/9)^2-v(t), 
q(t) = 1/2*1/(((u(t)+400/9)^2)^(1/2))-1/20-w(t)/(1/2*((u(t)+400/9)^2)^(1/2)-1/40*(u(t)+400/9)^2-1/2*v(t))
];

I don't know how to find the Jordan Normal form of this system. But the phase diagram of the 2-D approximation of the transformed system is:

[SERVER UNABLE TO UPLOAD FILES, FIGURE WILL BE UPLOADED LATER BELOW]

My interpretation of the figure is that the transformed system is stable. Do you agree?

What should I believe, is there a one-dimensional manifold that leads to this system?

N.B. x(t) is a state variable, while c(t) is a control variable and q(t) an adjoint variable, so that initial values c(0) and q(0) can be chosen so as to guide x(t) to the stationary state, if possible.

My complete code is below:

> A:=1: d:=1/20: s:=1/2: a:=1/2:
> b1:=2/100: b2:=3/100: w2:=1/2: w1:=1-w2: b:=w1*b1+w2*b2: b:
> f := x-> A*x^(a)-d*x:
> fx := x-> A*a*x^(a-1)-d:
> g := c-> s*c:
# Differential Equation System
> xdot := diff(x(t),t) = f(x(t))-c(t);
> cdot := diff(c(t),t) = g(c(t))*(fx(x(t))-q(t));
> qdot := diff(q(t),t) = -(q(t)-b1)*(q(t)-b2)+(q(t)-b)*rhs(cdot)/rhs(xdot);
# Stationary State (critical values of interest)
> xd:=eval(rhs(xdot),{x(t)=x, c(t)=c, q(t)=q}):
> cd:=eval(rhs(cdot),{x(t)=x, c(t)=c, q(t)=q}):
> qd:=eval(rhs(qdot),{x(t)=x, c(t)=c, q(t)=q}):
> qs:=b: qs:=convert(%,rational):
> xs:=fsolve(eval(cd/c,q=qs),x): xs:=convert(%,rational):
> cs:=fsolve(eval(xd,x=xs),c): cs:=convert(%,rational):
> ss:=evalf([xs,cs,qs]): convert(%,rational);
>
> #TRANSFORMATION 1
> newsys:= [
> diff(x(t),t) = rhs(xdot)*rhs(xdot),
> diff(c(t),t) = rhs(cdot)*rhs(xdot),
> diff(q(t),t) = rhs(qdot)*rhs(xdot)
> ];
>
> newvars := [
> x(t) = u(t) + xs,
> c(t) = f(u(t)+xs) - v(t),
> q(t) = fx(u(t)+xs) - w(t) / g(f(u(t)+xs)-v(t))
> ];
> uvwsys := op(simplify(solve(eval(newsys,newvars),[diff(u(t),t), diff(v(t),t), diff(w(t),t)])) assuming u(t)>0);
> eval(uvwsys,{u=0,v=0,w=0});
> F := subs(u(t)=u,v(t)=v,w(t)=w,map(rhs,uvwsys)):
> VectorCalculus[Jacobian](F,[u,v,w]) assuming u>0:
> J := simplify(eval(VectorCalculus[Jacobian](F,[u,v,w]),{u=0,v=0,w=0}));
> LinearAlgebra:-JordanForm(J);
> F2 := eval(F[1..2],w=0);
> Fv:= expand(F2/v);
> Jv:= simplify(eval(VectorCalculus[Jacobian](Fv,[u,v]),[u=0,v=0]));
> evalf(LinearAlgebra[Eigenvectors](Jv));
> with(plots):
> fieldplot(F2,u=-.005 .. 0.005, v = -.0001 .. .0001, fieldstrength=fixed);
>
> #TRANSFORMATION 2
> newsys:= [
> diff(x(t),t) = rhs(xdot)*rhs(xdot),
> diff(c(t),t) = rhs(cdot)*rhs(xdot),
> diff(q(t),t) = rhs(qdot)*rhs(xdot)
> ];
>
> newvars := [
> x(t) = (u(t) + xs)^2,
> c(t) = f((u(t) + xs)^2) - v(t),
> q(t) = fx((u(t) + xs)^2) - w(t) / g(f((u(t) + xs)^2) - v(t))
> ];
> uvwsys := op(simplify(solve(eval(newsys,newvars),[diff(u(t),t), diff(v(t),t), diff(w(t),t)])) assuming u(t)>0);
> eval(uvwsys,{u=0,v=0,w=0});
> F := subs(u(t)=u,v(t)=v,w(t)=w,map(rhs,uvwsys)):
> VectorCalculus[Jacobian](F,[u,v,w]) assuming u>0:
> J := simplify(eval(VectorCalculus[Jacobian](F,[u,v,w]),{u=0,v=0,w=0}));
> LinearAlgebra:-JordanForm(J);
> F2 := eval(F[1..2],w=0);
> Fv:= expand(F2/v);
> Jv:= simplify(eval(VectorCalculus[Jacobian](Fv,[u,v]),[u=0,v=0]));
> evalf(LinearAlgebra[Eigenvectors](Jv));
> with(plots):
> fieldplot(F2,u=-.005 .. 0.005, v = -.0001 .. .0001, fieldstrength=fixed);
>
> #TRANSFORMATION 3
> newsys:= [
> diff(x(t),t) = rhs(xdot)*rhs(xdot)^2,
> diff(c(t),t) = rhs(cdot)*rhs(xdot)^2,
> diff(q(t),t) = rhs(qdot)*rhs(xdot)^2
> ];
>
> newvars := [
> x(t) = (u(t) + xs)^2,
> c(t) = f((u(t) + xs)^2) - v(t),
> q(t) = fx((u(t) + xs)^2) - w(t) / g(f((u(t) + xs)^2) - v(t))
> ];
> uvwsys := op(simplify(solve(eval(newsys,newvars),[diff(u(t),t), diff(v(t),t), diff(w(t),t)])) assuming u(t)>0);
> eval(uvwsys,{u=0,v=0,w=0});
> F := subs(u(t)=u,v(t)=v,w(t)=w,map(rhs,uvwsys)):
> VectorCalculus[Jacobian](F,[u,v,w]) assuming u>0:
> J := simplify(eval(VectorCalculus[Jacobian](F,[u,v,w]),{u=0,v=0,w=0}));
> LinearAlgebra:-JordanForm(J);
> F2 := eval(F[1..2],w=0);
> Fv:= expand(F2/v^2);
> Jv:= simplify(eval(VectorCalculus[Jacobian](Fv,[u,v]),[u=0,v=0]));
> evalf(LinearAlgebra[Eigenvectors](Jv));
> with(plots):
> fieldplot(F2,u=-.005 .. 0.005, v = -.0001 .. .0001, fieldstrength=fixed);
>
 

As I noted, a typo crept up in the system analyzed by Robert. Two days ago I reported my inability to diagonalize the system once the typo had been corrected. I have now found a transformation that diagonalizes the system. I report on it below. I'm still unsure of the interpretation of it all though...

Using a transformation similar to the one Robert applied to the typo-contaminated system (TRANSFORMATION 1 below) yields the non-diagonal matrix:

But the following transformation works:

newsys := [
diff(x(t),t) = (x(t)^(1/2)-1/20*x(t)-c(t))^2, 
diff(c(t),t) = 1/2*c(t)*(1/2*1/(x(t)^(1/2))-1/20-q(t))*(x(t)^(1/2)-1/20*x(t)-c(t)), 
diff(q(t),t) = (-(q(t)-1/50)*(q(t)-3/100)
+1/2*(q(t)-1/40)*c(t)*(1/2*1/(x(t)^(1/2))-1/20-q(t))/(x(t)^(1/2)-1/20*x(t)-c(t)))*(x(t)^(1/2)-1/20*x(t)-c(t))
];

newvars := [
x(t) = (u(t)+400/9)^2, 
c(t) = ((u(t)+400/9)^2)^(1/2)-1/20*(u(t)+400/9)^2-v(t), 
q(t) = 1/2*1/(((u(t)+400/9)^2)^(1/2))-1/20-w(t)/(1/2*((u(t)+400/9)^2)^(1/2)-1/40*(u(t)+400/9)^2-1/2*v(t))
];

Its Jordan form is diagonal:

MATRIX([[0, 0, 0], [0, -187/108, 0], [0, 0, 0]]);

The phase diagram of the 2-D approximation of the transformed system is:

[SERVER UNABLE TO UPLOAD FILES, FIGURE WILL BE UPLOADED LATER BELOW]

My interpretation of the phase diagram is that the stationary state (u=v=w=0) is unstable. Do you agree?

But now consider a different transformation, which allows the derivative of x(t) to be both positive and negative,

newsys := [
diff(x(t),t) = (x(t)^(1/2)-1/20*x(t)-c(t))^3, 
diff(c(t),t) = 1/2*c(t)*(1/2*1/(x(t)^(1/2))-1/20-q(t))*(x(t)^(1/2)-1/20*x(t)-c(t))^2, 
diff(q(t),t) = (-(q(t)-1/50)*(q(t)-3/100)
+1/2*(q(t)-1/40)*c(t)*(1/2*1/(x(t)^(1/2))-1/20-q(t))/(x(t)^(1/2)-1/20*x(t)-c(t)))*(x(t)^(1/2)-1/20*x(t)-c(t))^2
];

newvars := [
x(t) = (u(t)+400/9)^2, 
c(t) = ((u(t)+400/9)^2)^(1/2)-1/20*(u(t)+400/9)^2-v(t), 
q(t) = 1/2*1/(((u(t)+400/9)^2)^(1/2))-1/20-w(t)/(1/2*((u(t)+400/9)^2)^(1/2)-1/40*(u(t)+400/9)^2-1/2*v(t))
];

I don't know how to find the Jordan Normal form of this system. But the phase diagram of the 2-D approximation of the transformed system is:

[SERVER UNABLE TO UPLOAD FILES, FIGURE WILL BE UPLOADED LATER BELOW]

My interpretation of the figure is that the transformed system is stable. Do you agree?

What should I believe, is there a one-dimensional manifold that leads to this system?

N.B. x(t) is a state variable, while c(t) is a control variable and q(t) an adjoint variable, so that initial values c(0) and q(0) can be chosen so as to guide x(t) to the stationary state, if possible.

My complete code is below:

> A:=1: d:=1/20: s:=1/2: a:=1/2:
> b1:=2/100: b2:=3/100: w2:=1/2: w1:=1-w2: b:=w1*b1+w2*b2: b:
> f := x-> A*x^(a)-d*x:
> fx := x-> A*a*x^(a-1)-d:
> g := c-> s*c:
# Differential Equation System
> xdot := diff(x(t),t) = f(x(t))-c(t);
> cdot := diff(c(t),t) = g(c(t))*(fx(x(t))-q(t));
> qdot := diff(q(t),t) = -(q(t)-b1)*(q(t)-b2)+(q(t)-b)*rhs(cdot)/rhs(xdot);
# Stationary State (critical values of interest)
> xd:=eval(rhs(xdot),{x(t)=x, c(t)=c, q(t)=q}):
> cd:=eval(rhs(cdot),{x(t)=x, c(t)=c, q(t)=q}):
> qd:=eval(rhs(qdot),{x(t)=x, c(t)=c, q(t)=q}):
> qs:=b: qs:=convert(%,rational):
> xs:=fsolve(eval(cd/c,q=qs),x): xs:=convert(%,rational):
> cs:=fsolve(eval(xd,x=xs),c): cs:=convert(%,rational):
> ss:=evalf([xs,cs,qs]): convert(%,rational);
>
> #TRANSFORMATION 1
> newsys:= [
> diff(x(t),t) = rhs(xdot)*rhs(xdot),
> diff(c(t),t) = rhs(cdot)*rhs(xdot),
> diff(q(t),t) = rhs(qdot)*rhs(xdot)
> ];
>
> newvars := [
> x(t) = u(t) + xs,
> c(t) = f(u(t)+xs) - v(t),
> q(t) = fx(u(t)+xs) - w(t) / g(f(u(t)+xs)-v(t))
> ];
> uvwsys := op(simplify(solve(eval(newsys,newvars),[diff(u(t),t), diff(v(t),t), diff(w(t),t)])) assuming u(t)>0);
> eval(uvwsys,{u=0,v=0,w=0});
> F := subs(u(t)=u,v(t)=v,w(t)=w,map(rhs,uvwsys)):
> VectorCalculus[Jacobian](F,[u,v,w]) assuming u>0:
> J := simplify(eval(VectorCalculus[Jacobian](F,[u,v,w]),{u=0,v=0,w=0}));
> LinearAlgebra:-JordanForm(J);
> F2 := eval(F[1..2],w=0);
> Fv:= expand(F2/v);
> Jv:= simplify(eval(VectorCalculus[Jacobian](Fv,[u,v]),[u=0,v=0]));
> evalf(LinearAlgebra[Eigenvectors](Jv));
> with(plots):
> fieldplot(F2,u=-.005 .. 0.005, v = -.0001 .. .0001, fieldstrength=fixed);
>
> #TRANSFORMATION 2
> newsys:= [
> diff(x(t),t) = rhs(xdot)*rhs(xdot),
> diff(c(t),t) = rhs(cdot)*rhs(xdot),
> diff(q(t),t) = rhs(qdot)*rhs(xdot)
> ];
>
> newvars := [
> x(t) = (u(t) + xs)^2,
> c(t) = f((u(t) + xs)^2) - v(t),
> q(t) = fx((u(t) + xs)^2) - w(t) / g(f((u(t) + xs)^2) - v(t))
> ];
> uvwsys := op(simplify(solve(eval(newsys,newvars),[diff(u(t),t), diff(v(t),t), diff(w(t),t)])) assuming u(t)>0);
> eval(uvwsys,{u=0,v=0,w=0});
> F := subs(u(t)=u,v(t)=v,w(t)=w,map(rhs,uvwsys)):
> VectorCalculus[Jacobian](F,[u,v,w]) assuming u>0:
> J := simplify(eval(VectorCalculus[Jacobian](F,[u,v,w]),{u=0,v=0,w=0}));
> LinearAlgebra:-JordanForm(J);
> F2 := eval(F[1..2],w=0);
> Fv:= expand(F2/v);
> Jv:= simplify(eval(VectorCalculus[Jacobian](Fv,[u,v]),[u=0,v=0]));
> evalf(LinearAlgebra[Eigenvectors](Jv));
> with(plots):
> fieldplot(F2,u=-.005 .. 0.005, v = -.0001 .. .0001, fieldstrength=fixed);
>
> #TRANSFORMATION 3
> newsys:= [
> diff(x(t),t) = rhs(xdot)*rhs(xdot)^2,
> diff(c(t),t) = rhs(cdot)*rhs(xdot)^2,
> diff(q(t),t) = rhs(qdot)*rhs(xdot)^2
> ];
>
> newvars := [
> x(t) = (u(t) + xs)^2,
> c(t) = f((u(t) + xs)^2) - v(t),
> q(t) = fx((u(t) + xs)^2) - w(t) / g(f((u(t) + xs)^2) - v(t))
> ];
> uvwsys := op(simplify(solve(eval(newsys,newvars),[diff(u(t),t), diff(v(t),t), diff(w(t),t)])) assuming u(t)>0);
> eval(uvwsys,{u=0,v=0,w=0});
> F := subs(u(t)=u,v(t)=v,w(t)=w,map(rhs,uvwsys)):
> VectorCalculus[Jacobian](F,[u,v,w]) assuming u>0:
> J := simplify(eval(VectorCalculus[Jacobian](F,[u,v,w]),{u=0,v=0,w=0}));
> LinearAlgebra:-JordanForm(J);
> F2 := eval(F[1..2],w=0);
> Fv:= expand(F2/v^2);
> Jv:= simplify(eval(VectorCalculus[Jacobian](Fv,[u,v]),[u=0,v=0]));
> evalf(LinearAlgebra[Eigenvectors](Jv));
> with(plots):
> fieldplot(F2,u=-.005 .. 0.005, v = -.0001 .. .0001, fieldstrength=fixed);
>
 

sorry about all the white space above, here's the system again :

newsys:=[
diff(x(t),t) = (x(t)^(1/2)-c(t))^3, 
diff(c(t),t) = c(t)*(1/2*1/(x(t)^(1/2))-q(t))*(x(t)^(1/2)-c(t))^2, 
diff(q(t),t) = (-(q(t)-1/50)*(q(t)-3/100)
+(q(t)-1/40)*c(t)*(1/2*1/(x(t)^(1/2))-q(t))
/(x(t)^(1/2)-c(t)))*(x(t)^(1/2)-c(t))^2
];

and the suggested transformation:

newvars:={
c(t) = -v(t)+(u(t)+400)^(1/2), 
q(t) = -w(t)/(-v(t)+(u(t)+400)^(1/2))+1/2*1/((u(t)+400)^(1/2)), 
x(t) = u(t)+400
};

sorry about all the white space above, here's the system again :

newsys:=[
diff(x(t),t) = (x(t)^(1/2)-c(t))^3, 
diff(c(t),t) = c(t)*(1/2*1/(x(t)^(1/2))-q(t))*(x(t)^(1/2)-c(t))^2, 
diff(q(t),t) = (-(q(t)-1/50)*(q(t)-3/100)
+(q(t)-1/40)*c(t)*(1/2*1/(x(t)^(1/2))-q(t))
/(x(t)^(1/2)-c(t)))*(x(t)^(1/2)-c(t))^2
];

and the suggested transformation:

newvars:={
c(t) = -v(t)+(u(t)+400)^(1/2), 
q(t) = -w(t)/(-v(t)+(u(t)+400)^(1/2))+1/2*1/((u(t)+400)^(1/2)), 
x(t) = u(t)+400
};

The change of speed suggested by Robert appears to be problematic,

with the new referential, the derivative of x(t) with respect to the new time index is now always positive, but I would like to understand the behavior of the system for both positive and negative derivatives (speeds) of x(t). I'm not sure how to do that when the right-hand side of x'(t) is always positive,

for that reason, I had suggested another change of speed in my earlier post on scimath (see reference below), which yields a cube rather than a square on the right hand side, allowing us to study both positive and negative derivatives.

the reason for this preferred change of speed is that now the old and the new time scales always move in the same direction -- with the change of speed suggested by Robert, they move in the same direction if x is rising, but in opposite directions if x is falling.

unfortunately I'm still unable to diagonalize or reduce the system (whether with squares or with cubes). As I explained earlier, there was a typo in the system analyzed by Robert (a missing square root) -- once the typo is corrected the suggested method for diagonalization fails. I'm stuck here. I fear the system may not be diagonalized.

This is, I think, the system I want to analyze:

newsys:= [
> diff(x(t),t) = rhs(xdot)*rhs(xdot)^2,
> diff(c(t),t) = rhs(cdot)*rhs(xdot)^2,
> diff(q(t),t) = rhs(qdot)*rhs(xdot)^2
> ];

Spurred on by Robert's earlier advice, I've transformed the system by setting v=x'(t) and w=c'(t) and analyzed (x,v,w), but I don't know how to diagonalize the system.

I think that this system is characterized by a center manifold and the equilibrium can be approached only from one side (for x rising), but this is based on intuition (detailed below) rather than a rigorous analysis, so any hints on how to proceed from now will be greatly appreciated.

The change of speed suggested by Robert appears to be problematic,

with the new referential, the derivative of x(t) with respect to the new time index is now always positive, but I would like to understand the behavior of the system for both positive and negative derivatives (speeds) of x(t). I'm not sure how to do that when the right-hand side of x'(t) is always positive,

for that reason, I had suggested another change of speed in my earlier post on scimath (see reference below), which yields a cube rather than a square on the right hand side, allowing us to study both positive and negative derivatives.

the reason for this preferred change of speed is that now the old and the new time scales always move in the same direction -- with the change of speed suggested by Robert, they move in the same direction if x is rising, but in opposite directions if x is falling.

unfortunately I'm still unable to diagonalize or reduce the system (whether with squares or with cubes). As I explained earlier, there was a typo in the system analyzed by Robert (a missing square root) -- once the typo is corrected the suggested method for diagonalization fails. I'm stuck here. I fear the system may not be diagonalized.

This is, I think, the system I want to analyze:

newsys:= [
> diff(x(t),t) = rhs(xdot)*rhs(xdot)^2,
> diff(c(t),t) = rhs(cdot)*rhs(xdot)^2,
> diff(q(t),t) = rhs(qdot)*rhs(xdot)^2
> ];

Spurred on by Robert's earlier advice, I've transformed the system by setting v=x'(t) and w=c'(t) and analyzed (x,v,w), but I don't know how to diagonalize the system.

I think that this system is characterized by a center manifold and the equilibrium can be approached only from one side (for x rising), but this is based on intuition (detailed below) rather than a rigorous analysis, so any hints on how to proceed from now will be greatly appreciated.

restart: with(plots):

f:=x->m*log(x)/2^m+(1-x^m)/(1+x)^m;

                                               m
                              m log(x)    1 - x
                    f := x -> -------- + --------
                                  m             m
                                 2       (1 + x)

> fx:=(x,m)->D(f)(x);

                       fx := (x, m) -> D(f)(x)

> plot3d(fx(x,m),x=1..10,m=0..10, axes=boxed);

And this may provide further clues:

fxm:=(x,m)->diff(fx(x,m),m);

plot3d(fxm(x,m),x=1..10,m=0..10, axes=boxed);

restart: with(plots):

f:=x->m*log(x)/2^m+(1-x^m)/(1+x)^m;

                                               m
                              m log(x)    1 - x
                    f := x -> -------- + --------
                                  m             m
                                 2       (1 + x)

> fx:=(x,m)->D(f)(x);

                       fx := (x, m) -> D(f)(x)

> plot3d(fx(x,m),x=1..10,m=0..10, axes=boxed);

And this may provide further clues:

fxm:=(x,m)->diff(fx(x,m),m);

plot3d(fxm(x,m),x=1..10,m=0..10, axes=boxed);