Hint: What's equivalent to the remainder of p(x)/(x-a) for a polynomial p?

@Thomas Dean As I said above, it's documented at ?evalapply.

@Thomas Dean Yes, [f,g](a,b) becomes [f(a,b), g(a,b)]. It's as simple as that. In terms of a theory of operators, the function application operator distributes over the list-building operator.

I think that he means the syntax by which [cos,sin](t) becomes [cos(t),sin(t)]. If that's what you want, it's documented at ?evalapply, but there isn't much to say about it: It works that way, and that's it.

Another possibility is the prefix-form list constructor `[]`. 

@goli  The command textplot3d does not support any directional or rotational options, even in Maple 2020.

However, the colored-dots Answer below can be easily adapted to 3D plots, even in Maple 15.

@Carl Love I just posted a new Answer that doesn't use textplot.

@goli Sorry, there's no rotation option for textplot in Maple 15. This is a situation where you'll need a newer version.

I can concoct a solution that uses colored symbols (that are not arrowheads) to indicate the direction of travel.

@goli Did you re-execute the entire worksheet by clicking on the !!! on the toolbar?

@hym100 So, my solution is

C:= solve(int(c*x^2, x= 0..3)=1);
int(piecewise(x < 0, 0, x <= 3, C*x^2, 0), x= -1..1);
 

@goli Change the procedure Ang to

Ang:= tau->
    eval(
        eval(
            arctan(diff(y(t),t)*(Zmax-Zmin), diff(z(t),t)*(Ymax-Ymin)),
            sys11
        ),
        dsol(tau)
    )
:

I also recommend that you remove the stepsize option from DEplot and change the t-range to -3..3.

@goli I'll rewrite it to use plain eval. It should only take a few minutes.

@goli Check your Maple help (?eval) to see if there's an option eval[recurse]. If there is, it'll be mentioned right at the top of the help page.

@Ramakrishnan Yes, it's a great example: It shows using dsolve with parameters, sophisticated use of eval, and two-argument arctan.

But the OP, goli, was mistaken in calling these PDEs. They're actually ODEs. The distinction is the number of independent variables with respect to which derivatives are taken. In this case, there's only one, t. One such variable makes them ODEs; more than one, PDEs.

@goli Beginning exactly where your worksheet leaves off, do this:

PrePlot:= %:
(Zmin,Zmax):= (-1,1): (Ymin,Ymax):= (0,1):
dsol:= dsolve(
    {sys11[], z(0)=z0, y(0)=y0}, 
    numeric, parameters= [z0,y0]
):
Ang:= tau->
    eval[recurse](
        arctan(diff(y(t),t)*(Zmax-Zmin), diff(z(t),t)*(Ymax-Ymin)),
        [dsol(tau)[], sys11[]]
    )
:
Arrows:= proc(ics)
local k;
    dsol(parameters= eval([z(0),y(0)], ics));
    plots:-textplot(
        [seq](
            eval([z(t), y(t), ">", rotation= Ang(k)], dsol(k)),
            k= [-1, 0, 1]
        ),
        font= [Symbol, Bold, 30], _rest
    )
end proc
:
plots:-display(
    [PrePlot, Arrows(ics1, color= red), Arrows(ics2, color= blue)]
);

@MapleUser09 My procedure is not intended to compute Fibonacci numbers! I posted it because it uses a computational technique similar to what you'd use for Fibonacci numbers. You need to write that. The defining property of Fibonacci numbers is

Fib(1)=1, Fib(2)=1, F(n)=Fib(n-1)+Fib(n-2)