@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)

@MapleUser09 My code above is not for Fibonacci numbers. Once you've written your own Fib with option remember, you could say Fib(100), Fib(1000), Fib(10000), and all will be computed very quickly. I'm trying to help you without completely giving away the answers. These are very easy homework problems.

Couldn't you at least try to rewrite your homework problems in your own words?

What have you done? What's your Fib procedure?

@Ramakrishnan Add the option font= [Symbol, Bold, n] where n is some positive integer proportional to the desired size. But note that this number is not directly related to the thickness.

plots:-display(
    PrePlot,
    plots:-textplot(
        [seq([k, f(k), ">", rotation= Angle(k)], k= ceil(Xmin)..floor(Xmax))],
        color= red, font= [Symbol, Bold, 40]
    ),
    thickness= 5
);

Couldn't you at least try to rewrite your homework problems in your own words? and come up with a better title? Your style is lame.

I find that rewriting a problem often makes me realize the solution.

@Ramakrishnan Thank you so much for catching that. I've now corrected the Answer. As you've probably guessed, I originally used A in both places, but at the last moment changed the name in the definition to Angle for clarity.

@dharr The computational effort required for gridrefine=n is an exponential function of n in both time and memory. The highest that I've ever used was gridrefine=5.

@acer I get your point, and combined with the OPs later comments, I see that maximally consistent subsets aren't sufficient. 

Suppose that the equations were {a=1, a=2, b=3}. Then I think now that the desired subset is {b=3}. But I don't see how to generalize that to cases where the "good" variables and "bad" variables are mixed in the same equations. 

@emendes Would it be sufficient to find a maximally consistent subset?

You wrote, "otherwise maple needlessly spend time for finding the gcd." 

I don't believe that unless you can show me an example.

Given: polynomials p, q, and f where f is a previously known common divisor of p and q.

Algorithm A: Compute GCD(p/f, q/f).
Algorithm B: Compute GCD(p,q)/f.

Can you show me an example where A is faster than B?

@Ugurgozutok If the polynomial is not integer, then the situation may be hopeless. What kind of coefficients do you have? For decimal coefficients, you may be able to find an approximate factorization with the SNAP package (see ?SNAP). I don't have any experience using it with parameters (such as your d).

My algorithm above does not care whether the polynomial is linear in d, but it does require that both d and the coefficients be integer.

@Carl Love I just updated the example in my Answer. It's a 20-digit (67-bit) value of d found in 2.3 seconds.