Download DirGraph.mw

Why do you think that the desired solution must be obtained for _C1 = 0?
Use e.g.
S:=dsolve({E1, phi(0) = b*a/(a*c-b^2)},  phi(xi));
symplify([S], symbolic);

You will find your solution, but be careful about branches.

(Note also that the ODE is not in an explicite form, so even with an initial condition it could have infinitely many solutions).

Note first that A^n can be computed symbolically using LinearAlgebra:-MatrixPower(A,n);
The entries are the Fibonacci numbers, and as Rouben said, it is easy to see that the maximum entry is fibonacci(n+1) at (1,1).
Here is a procedure to compute the inverse Fibonacci numbers.

inversefib:=proc(y::posint)
local n1,n2, s1,s2;
if y=1 then return 'Exact'=2 fi;
n1:=ln( (1/2)*sqrt(5)*y+(1/2)*sqrt(5*y^2+4) )/ln((sqrt(5)+1)/2); # n even
n2:=ln( (1/2)*sqrt(5)*y+(1/2)*sqrt(5*y^2-4) )/ln((sqrt(5)+1)/2); # n odd
s1,s2:=simplify~([n1,n2])[];
if   s1::even then  'Exact'=s1
elif s2::odd  then  'Exact'=s2
else  'Next'= ceil(s1) fi
end proc:

Examples.
inversefib(2019);
     Next = 18
inversefib(3524578);
     Exact = 33
inversefib(10^12);
     Next = 60

remove[flatten](member, A, B);

If you want it inplace, use

remove[flatten,inplace](member, A, B);

plot(arctan(x),x=-20..20, 
title=typeset(
H__0,": ", sigma[1]^2=sigma[2]^2, "     ",  
H__1,": ", sigma[1]^2<>sigma[2]^2, "\n", MyTitle)
);

Yes, it seems that

seq(C[i+1]-C[i], i = 1 .. N-1):
does not finish the job.
But,
seq(C[i+1]-i^2, i = 1 .. N-1):
works.
Probably a bug.

The limit does exist in C.
For a real function it is possible to use

RealDomain:-limit(sqrt(x^2 - 3 *x + 2), x=2, left);
        undefined

Or, execute first
with(RealDomain);

However I'd recommend RealDomain only if the user really does not know (yet) about complex numbers & functions.
(You should know that almost all mathematical functions are defined in C.)

FC:= proc(n) 
uses plots,plottools;
local a,b,u;
{seq(seq(a/b, a = 0 .. b), b = 1 .. n)};
display(seq(circle([u,1/2/denom(u)^2],1/2/denom(u)^2), u=%)
        ,scaling = constrained, color = blue, size = [800, 800])
end proc:

FC(6);

See also Apollonian circles

print~([a, b]);  means  [print(a), print(b)];

Now, print(a)  prints a (as a side effect) and returns NULL.
The last output is [NULL,NULL]. If you want to eliminate it, use print~([a, b]):  (colon instead of semicolon).

Digits:=150:
A:=Matrix(eval(M,P=0),shape=symmetric):
B:=Matrix(-diff(M,P), shape=symmetric):
E:=LinearAlgebra:-Eigenvalues(A,B):
E1:=Re(E);
E2:=Im(E); # =0
Esorted:=sort(E1);
evalf[15](Esorted[1..15]), evalf[15](Esorted[-15..-1]);

Yes, shape=triangular[upper, unit]  seems to have a serious problem.

The best (partial) workaround is probably

U:= Matrix(Matrix(4,4, shape=triangular[upper, unit], symbol =u), shape=triangular[upper]);

(partial, because the diagonal elements are now assignable.)

s:=solve((hv*hw+lv*lw)/(lv+hv)*(lv+hv) = dw*(lv+hv), hv):
sign(s)*numer(s)/sort(expand((sign(s)*denom(s))),order=plex(hw,dw));
#   Do you think it's worth the effort?

solve({seq(a*v + b*u - w)}, {a,b});

(If w is not in the span of u,v, the result is NULL).