Question: approximating a center manifold

I have a 3-dimensional dynamic system of ODEs characterized by zero eigenvalues at the origin. I would like to know if it's possible to find one trajectory converging towards the origin.

My system is derived from a control problem and has consequently control-like and state-like variables -- while the neighborhood is typically unstable I need to know if an appropriate choice of initial values can steer the system towards the origin;  I expect at most one such trajectory.

I'm trying to obtain a center manifold reduction near the origin of my system, to establish whether the center manifold is stable or unstable. Any help would be greatly appreciated.

Here's the little I know about center manifold theory: If a dynamic system has zero eigenvalue(s) near a critical point of interest, the Jacobian is uninformative about the local dynamics. To apply the center manifold theorem, one first transforms the system using the eigenvectors. For instance, if V denotes the matrix of eigenvectors, I should transform my [u,v,w] system into [uu,vv,ww] as follows:

[u,v,w] = V. [uu,vv,ww];

I get stuck at this point because my calculations yield that one of the eigenvectors is (0,0,0). I'm not sure how to proceed from there.

My system and my analysis is given below. The last line shows the calculations that yield the eigenvectors.

Note:  The system named sys is the system of interest, with the system newsys obtained by a change of speed. I'd like to apply a center manifold reduction to newsys. The variables are [u,v,w] and the critical point of interest is the origin [0,0,0]; eigenvalues are zero at the origin.

restart;
# Parameter and Functional Assignments
params := [ A=2, d=1/20, s=1, a=1/2, b1=2/100, b2=3/100, 
w1=1/2, w2=1/2, b=1/40 ];

f   := x -> A*x^(a)-d*x:
fx  := x -> A*a*x^(a-1)-d:
fxx := x -> A*a*(a-1)*x^(a-2):
ifx := y -> (A*a/(d+y))^(1-a):
h   := y -> (fxx(ifx(y)))*(f(ifx(y))): h(y):

# Differential Equation System
udot := diff(u(t),t) = (1-u(t))*(u(t)*(b+v(t))-s*(v(t)+w(t))):
vdot := diff(v(t),t) = u(t)*h(b+v(t)):
wdot := diff(w(t),t) = -(w(t)+b1-b)*(w(t)+b2-b)+s*w(t)*(v(t)
+w(t))/u(t)*(1-u(t)):
sys:=[udot,vdot,wdot]: eval(sys,params);

newsys := subs(u(t)=u,v(t)=v,w(t)=w,map(rhs,sys *~ u(t))):
newsys := simplify(newsys) assuming positive:
eval(newsys,params);

mtaylor(newsys[2],[u,v,w],3):
newsystaylor := [expand(newsys[1]),%,newsys[3]]:
eval(newsystaylor,params);

Jac := VectorCalculus:-Jacobian(newsystaylor,[u,v,w]):
Jac := eval(simplify(Jac),{u*v=0,u*w=0,v*w=0,u^2=0,v^2=0,w^2=0,u^3=0}): 
eval(Jac,params);

# For an Order-3 Approximation
eval(Jac[2,1]/u,[u=0,v=0,w=0]):
C3 := -2*A*a*(a-1)*(A^(-2*a^2+4*a-1)*a^(-2*(a-1)^2)
*(d^2+2*d*b+b^2)^(a*(a-2))*d^2-A^(-(a-1)^2)*a^(-(a-1)^2)
*d*(d^2+2*d*b+b^2)^(-a)*(b+d)^(-1+a^2)*b^2-A^(-(a-1)^2)
*a^(-(a-1)^2)*d^3*(d^2+2*d*b+b^2)^(-a)*(b+d)^(-1+a^2)
+2*A^(-2*a^2+4*a-1)*a^(-2*(a-1)^2)*(d^2+2*d*b+b^2)^(a*(a-2))
*d*b+A^(-2*a^2+4*a-1)*a^(-2*(a-1)^2)*(d^2+2*d*b+b^2)^(a*(a-2))
*b^2-2*A^(-(a-1)^2)*a^(-(a-1)^2)*d^2*(d^2+2*d*b+b^2)^(-a)*(b+d)^(-1+a^2)*b):
 
is(Jac[2,1]=-C3*u);

Jac[2,1] := -C*u: Jac;

e := LinearAlgebra:-Eigenvalues(Jac): 
eval(e[1],[u=0,v=0]);
eval(e[2],[u=0,v=0]);
eval(e[3],[u=0,v=0]); 

coulditbe( e[2] in 'SetOf'('real') );

Jac0 := eval(VectorCalculus:-Jacobian(newsystaylor,[u,v,w]),{u=0,v=0,w=0}):

C2 := -(b-b1)*(b-b2): expand(%):
is(Jac0[3,1]=C2):
Jac0[3,1] := B: Jac0;

(E, V) := LinearAlgebra:-Eigenvectors(Jac0);
(J, Q) := LinearAlgebra:-JordanForm(Jac0,output=['J','Q'] );

                            [0    0    0]
                            [           ]
                            [0    0    0]
                            [           ]
                            [B    0    0]


                              [0]  [0    0    0]
                              [ ]  [           ]
                      E, V := [0], [0    1    0]
                              [ ]  [           ]
                              [0]  [1    0    0]


                         [0    1    0]  [0    1    0]
                         [           ]  [           ]
                 J, Q := [0    0    0], [0    1    1]
                         [           ]  [           ]
                         [0    0    0]  [B    0    0]