But applied is documented: ?type,applied.  What's not documented is the reason for it.  It's not mentioned in the ?updates,Maple9_5 pages, where new features should be introduced.

Where did you find ComplementaryIndices?  It seems to be undocumented in the Maple help system.

Where did you find ComplementaryIndices?  It seems to be undocumented in the Maple help system.

It's really unfortunate that the Maple type function refers to what should really be a "function call".  Since this goes way back to the very beginnings of Maple, I guess we're stuck with it to preserve backward compatibility.  The type applied  is identical to function, however, so in principle I think all references to function could be replaced by applied.  I don't know why applied was introduced (some time after Maple 6), if not for this purpose.

It's really unfortunate that the Maple type function refers to what should really be a "function call".  Since this goes way back to the very beginnings of Maple, I guess we're stuck with it to preserve backward compatibility.  The type applied  is identical to function, however, so in principle I think all references to function could be replaced by applied.  I don't know why applied was introduced (some time after Maple 6), if not for this purpose.

Click on the green up-arrow key, then follow the instructions in the File Manager window that should pop up.

If you use one of the classical methods with two step sizes, you can use Richardson extrapolation to try to estimate the error. 

If you use one of the classical methods with two step sizes, you can use Richardson extrapolation to try to estimate the error. 

One problem seems to be that H is an expression, not a function, so H(u,v) doesn't make any sense.  So it would have to be color=H, not color=H(u,v).  But there's an additional problem: H also contains ν_x and ν_z., which plot3d knows nothing about.

Since I don't know much about differential geometry or the DifferentialGeometry and Tensor packages, I don't think I can help with the rest.

One problem seems to be that H is an expression, not a function, so H(u,v) doesn't make any sense.  So it would have to be color=H, not color=H(u,v).  But there's an additional problem: H also contains ν_x and ν_z., which plot3d knows nothing about.

Since I don't know much about differential geometry or the DifferentialGeometry and Tensor packages, I don't think I can help with the rest.

Please write exactly what you mean.  First you wrote f(g), now it's f(x)-g(x), but this does not really fit with the rest of what you wrote.  What do you want on which axis?

Please write exactly what you mean.  First you wrote f(g), now it's f(x)-g(x), but this does not really fit with the rest of what you wrote.  What do you want on which axis?

OK, I think this is the situation. 

The system, with the change of speed I defined in my previous posting, was

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

We are looking for solutions with x(t) -> 400/9, Q(t) -> 1/40, C(t) -> 40/9 as t -> 0 (but x(t) not constant).
I think there should be two such solutions on each side of x(t) = 400/9.

Change variables again, as follows.

> newvars:= {x(t)=u(t)+400/9,
        Q(t)=w(t)-3/17000*v(t)+1/40,
        C(t)=v(t) + sqrt(u(t)+400/9) - 1/20*(u(t) + 400/9)};
  uvwsys:= op(simplify(solve(eval(newsys,newvars),[diff(u(t),t), diff(v(t),t), diff(w(t),t)])));
  F:= subs(u(t)=u,v(t)=v,w(t)=w,map(rhs,uvwsys));

In these variables we want u(t), v(t), w(t) all to go to 0 (but u(t) not constant).  The factor 3/17000 was chosen so that the linearization of this system at the critical point (u,v,w) = (0,0,0) is diagonal.

> J:= simplify(eval(VectorCalculus[Jacobian](F,[u,v,w]),{u=0,v=0,w=0}));

                                   [0  0    0]
                                   [         ]
                                   [0  0    0]
                              J := [         ]
                                   [      -17]
                                   [0  0  ---]
                                   [      120]

The consequence of this is that as t -> +infinity, solutions approaching the critical point should have w going rapidly to 0, while u and v go to 0 on a slower time scale.  The behaviour should, I think, be like that of the system for u(t) and v(t) with w(t) = 0.  On the other hand, for t -> -infinity, most solutions will have w going rapidly away from 0, either to the positive or the negative side, but between these there should be a manifold of solutions with w -> 0.  On that manifold, again the behaviour should (I hope) be like that of the system for u(t) and v(t) with w(t) = 0.  Here is that system.
 

> F2:= eval(F[1..2],w=0);

F2 := [v^2, v*(-816000000+918000*u*sqrt(9*u+400)-1620*u*sqrt(9*u+400)*v+40800000*sqrt(9*u+400)-13842000*v*sqrt(9*u+400)-22950000*u+4841100*u*v-4860*v^2*u-103275*u^2+243*v*u^2+214680000*v-216000*v^2)/(6120000*(9*u+400))]

If you divide by v, you get a field with a saddle point at (0,0).

> Fv:= expand(F2/v);
  Jv:= simplify(eval(VectorCalculus[Jacobian](Fv,[u,v]),[u=0,v=0]));

                                   [   0      1]
                                   [           ]
                             Jv := [ 3    -259 ]
                                   [----  -----]
                                   [1600  10200]

>  evalf(LinearAlgebra[Eigenvectors](Jv));

                [ 0.03242808813]  [30.83746399  -17.29498033]
                [              ], [                         ]
                [-0.05782024499]  [         1.            1.]

However, note that for v < 0 this reverses the direction of the flow. For the actual field F2, the picture looks like this:

> with(plots):
  fieldplot(F2,u=-.005 .. 0.005, v = -.0001 .. .0001, fieldstrength=fixed);

The result is that there are two trajectories approaching (0,0) as t -> infinity on the left side of the picture, and two approaching (0,0) as t -> -infinity on the right side of the picture, in the directions of the eigenvectors of Jv.

OK, I think this is the situation. 

The system, with the change of speed I defined in my previous posting, was

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

We are looking for solutions with x(t) -> 400/9, Q(t) -> 1/40, C(t) -> 40/9 as t -> 0 (but x(t) not constant).
I think there should be two such solutions on each side of x(t) = 400/9.

Change variables again, as follows.

> newvars:= {x(t)=u(t)+400/9,
        Q(t)=w(t)-3/17000*v(t)+1/40,
        C(t)=v(t) + sqrt(u(t)+400/9) - 1/20*(u(t) + 400/9)};
  uvwsys:= op(simplify(solve(eval(newsys,newvars),[diff(u(t),t), diff(v(t),t), diff(w(t),t)])));
  F:= subs(u(t)=u,v(t)=v,w(t)=w,map(rhs,uvwsys));

In these variables we want u(t), v(t), w(t) all to go to 0 (but u(t) not constant).  The factor 3/17000 was chosen so that the linearization of this system at the critical point (u,v,w) = (0,0,0) is diagonal.

> J:= simplify(eval(VectorCalculus[Jacobian](F,[u,v,w]),{u=0,v=0,w=0}));

                                   [0  0    0]
                                   [         ]
                                   [0  0    0]
                              J := [         ]
                                   [      -17]
                                   [0  0  ---]
                                   [      120]

The consequence of this is that as t -> +infinity, solutions approaching the critical point should have w going rapidly to 0, while u and v go to 0 on a slower time scale.  The behaviour should, I think, be like that of the system for u(t) and v(t) with w(t) = 0.  On the other hand, for t -> -infinity, most solutions will have w going rapidly away from 0, either to the positive or the negative side, but between these there should be a manifold of solutions with w -> 0.  On that manifold, again the behaviour should (I hope) be like that of the system for u(t) and v(t) with w(t) = 0.  Here is that system.
 

> F2:= eval(F[1..2],w=0);

F2 := [v^2, v*(-816000000+918000*u*sqrt(9*u+400)-1620*u*sqrt(9*u+400)*v+40800000*sqrt(9*u+400)-13842000*v*sqrt(9*u+400)-22950000*u+4841100*u*v-4860*v^2*u-103275*u^2+243*v*u^2+214680000*v-216000*v^2)/(6120000*(9*u+400))]

If you divide by v, you get a field with a saddle point at (0,0).

> Fv:= expand(F2/v);
  Jv:= simplify(eval(VectorCalculus[Jacobian](Fv,[u,v]),[u=0,v=0]));

                                   [   0      1]
                                   [           ]
                             Jv := [ 3    -259 ]
                                   [----  -----]
                                   [1600  10200]

>  evalf(LinearAlgebra[Eigenvectors](Jv));

                [ 0.03242808813]  [30.83746399  -17.29498033]
                [              ], [                         ]
                [-0.05782024499]  [         1.            1.]

However, note that for v < 0 this reverses the direction of the flow. For the actual field F2, the picture looks like this:

> with(plots):
  fieldplot(F2,u=-.005 .. 0.005, v = -.0001 .. .0001, fieldstrength=fixed);

The result is that there are two trajectories approaching (0,0) as t -> infinity on the left side of the picture, and two approaching (0,0) as t -> -infinity on the right side of the picture, in the directions of the eigenvectors of Jv.

Integer linear programming is NP-hard.  Maple's Optimization package handles it using a branch-and-bound method, but this may take a very long time in some cases.  Integer nonlinear programming may be even harder.  For some of these problems, there may be no method much better than brute force.  If you have 50 or so variables, it may be unrealistic to expect to find the optimal solution.