You cannot simplify to 0 because it is not 0. But it is a constant.
You can check this in Maple, computing

diff(expr, x);   # 0

where expr is your expression. Do not forget to use u(x) and v(x)  instead of u, v.
 

Tunnell:=proc(N::posint)
local A, a,t, n, isq;
n:=mul(map(`^`~ @ `op`, map(t->[t[1],irem(t[2],2)], ifactors(N)[2])));
isq:=proc(n) local r:=isqrt(trunc(n)); `if`(r^2>n,r-1,r) end;

A:=proc(a,b,c)
  local x,y,z, r:=0;
  for z from 0 to isq(n/c) do
  for y from 0 to isq((n - c*z^2)/b) do
    x:=(n-c*z^2-b*y^2)/a;
    if issqr(x) then r:=r+`if`(z=0,1,2)*`if`(y=0,1,2)*`if`(x=0,1,2) fi
  od od: 
r 
end;

if    n::odd  and 2*A(1,2,32)=A(1,2,8) 
   or n::even and 2*A(2,8,64)=A(2,8,16)
then true else false fi
end:

select(Tunnell, [seq(1 .. 200)]);

[5, 6, 7, 13, 14, 15, 20, 21, 22, 23, 24, 28, 29, 30, 31, 34, 37, 38, 39, 41, 45, 46, 47, 52, 53, 54, 55, 56, 60, 61, 62, 63, 65, 69, 70, 71, 77, 78, 79, 80, 84, 85, 86, 87, 88, 92, 93, 94, 95, 96, 101, 102, 103, 109, 110, 111, 112, 116, 117, 118, 119, 120, 124, 125, 126, 127, 133, 134, 135, 136, 137, 138, 141, 142, 143, 145, 148, 149, 150, 151, 152, 154, 156, 157, 158, 159, 161, 164, 165, 166, 167, 173, 174, 175, 180, 181, 182, 183, 184, 188, 189, 190, 191, 194, 197, 198, 199]

Edit. The proc was optimized for speed.

A module is very convenient and has advantages, but it is not mandatory here.

with(Statistics):

Levy:= (m,s) -> Distribution(
                      PDF=(x->piecewise(x<0, 0, sqrt(s/2/Pi)*exp(-s/2/(x-m))/(x-m)^(3/2))), 
                      CDF= (x->erfc(sqrt(s/2/(x-m)))), 
                      RandomSample = (N->m +~ s /~ ( Quantile~(Normal(0, 1), Sample(Uniform(1/2, 1), N)) )^~2)
                    ):

X := RandomVariable(Levy(m,s)):
PDF(X, x);
Sample(X, 5);
Y := RandomVariable(Levy(1,2)):
PDF(Y, x);
Sample(Y, 6);

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]);