@Thomas Dean I executed the worksheet mw.mw from my earlier post in Maple 2021.2 (on Ubuntu).  It produced the inline graph as expected.  There must be something else in your setup that prevents the inline display.

You code works fine here.  Check your settings.

Download mw.mw

@Christian Wolinski When calling type() or indets(), as when calling any ordinary function, the arguments are evaluated first before passing to the function.  For instance,
indets(5^3, anything^Non(integer));
type(sqrt(4), radical);
return empty and false, respectively, because 5^3 gets evaluated to 125 and sqrt(4) gets evaluated to 2.  Similarly,
indets((-1)^(1/2), anything^Non(integer));
type(sqrt(-1), radical);
return empty and false, respectively, because (-1)^(1/2) and sqrt(-1) get evaluated to I.

Differences from an earlier version of Maple may be attributable to changes in the automatic simplification rules.  You may or may not like those changes but I wouldn't call then wrong or irrational.

@acer You are right.  What I posted under the heading "Answer to my own question" does not do justice.  My purpose there was to show the gist of an idea.  What we really want is an extension of Maple's D operator that may be applied to non-scalar-valued procs.  The following does exactly that.  I intend to add a few illustrative examples and upload it as a Post.

DD := proc(e::procedure)
	local idx, n, u, params, i;
	idx := op(procname);          # differentiation multi-index
	n := nops([op](1, eval(e)));  # n = number of e's parameters
	params := u[i] $i=1..n;
	if type(procname, indexed) then
		# all is well
	elif n = 1 then 
		idx := 1;
	else
		error "differentiate with resepct to which variable?";
	end;

	if type(e(params), scalar) then
		return D[idx](e);
	else
		unapply(map(x -> (D[idx](unapply(x, params)))(params), e(params)), params);
	end if;
end proc:

This proc extends Maple's D operator so that it applies to non-scalar-valued procedures.  For instance, if

z := (u,v) -> < a(u,v), b(u,v) >;

then DD[1,1,2](z) does what Maple's D[1,1,2](z) would do to a scalar-valued proc.
 

@Axel Vogt That's a good suggestion.  I am not as familiar with the VectorCalculus as I should be.  I will take a look.

I figured out the answer to my own question:

This generalizes nicely to arbitrary non-scalar-valued procedures.  I will explain the details in a Post.

@Carl Love That's illuminating, however I don't understand why you have And instead of Or in the defintion of Nparams.
Consider

f := proc(x, {y:=12})
	x -> x^2;
end proc:

Then Nparams(f) fails as written, but succeeds when we replace And with Or.

@acer that's excellent.  I had tried op(1,q) but that wansn't good.  It had not occurred to me to change q to eval(q).  Thanks for the lesson.

Regarding parameters vs arguments, I thought a bit about it when I posted the question.  To me it seems that the proper usages are "The procedure q has two parameters".  The procedure q takes two arguements" because arguments are passed to q.  But perhaps that's too much of a hair-splitting.

Hello @mmcdara, what you have proposed produces a vector whose first component is D[1](a(u,v)) which is not a very meaningful expression.  We want D[1](a)(u,v) instead.

@one man Getting rid of the denominators in the ODEs is not necessarily a good thing.  It's true that doing that results in neater looking equations, but that does not mean that it improves their solvability.

In the case of the current problem, after removing the denominators, the equations reduce to the deceptively pleasant-looking polynomial ODES:

diff(x(t),t) = z(t)*y(t)^2*(y(t)*z(t)+200)
diff(y(t),t) = -2*x(t)^2*z(t)*(x(t)*z(t)+100)
diff(z(t),t) = 20*(2*x(t)-y(t))*y(t)^2*x(t)^2

ICs := x(0) = 6.066066075, y(0) = -8.537156781, z(0) = 20  # as before

The solution of that system, however, behaves horribly.  The orbit still traces the same curve as before, as it should. However, having removed the denominators, the orbit is no longer a unit-speed curve.  The speed varies very sharply (almost discontinuously) in a couple of places, and as a result, the roots that we are seeking bunch up very close to each other.  So much so that fsolve() has difficulty in finding them.  Here is what the function f3 looks like along the orbit:

The 116 roots are hidden in that mess. Compare that to the corresponding graph in my original post.

By the way, in case you wish to try that yourself, the parameter range for a full orbit is t=0 .. 0.007604467601.

@mmcdara The purpose of my question is pedagogical.  In your solution you refer to L2 and Linfinity norms.  It seems to me that Otttimor does not know what those terms mean.  That''s why I asked him/her for an explanation.

@Otttimor Let's see if this helps.  Consider two lists of numbers:

vals1 := [1,3,5,6,6,6,7,5,3];
vals2 := [1,3,4,5,6,8,7,4,4];

Can you explain how you would want to measure the error/deviation between these lists? Your answer to that will determine the strategy to use in your problem.

By the way, you may plot that data through

dataplot([vals1,vals2], color=[red,blue]);

@Jonas111996 The coefficients of x1(t) and x2(t) on the right-hand sides of your differential equations are things like sin(t*sqrt(-3)) which are explicit functions of t.  That makes your system nonautonomous as Carl has already pointed out, so a phase portrait makes no sense.

Furthermore, due to the presence of the sqrt(-3), those coefficients are complex, making a 2-dimensional phase portrait impossible.

I suspect that you have gone off in the wrong direction in the process of arriving at your differential equations.  If you explain where those equations come from, someone should be able to orient you in the right direction.

You may want to ask the first author whose email is given in the article:

rxyao@sjtu.edu.cn

@acer I agree.  I wouldn't be surprised if a close look at the algorithm reveals better ways of achieving whatever that program is aiming to achieve.