The proof of the first formula is very simple but I rewrite F into F2 forreasons that will appear clearly during the the proof of your second formula.

Formula 1:

F := unapply(rsolve({F(1) = 1, F(2) = 1, F(n + 1) = F(n) + F(n - 1)}, F(n)), n):

rel2use := Phi = op([1, -1, 1], F(n)):
rel2use := { rel2use, C = op(1, F(n)) / rhs(rel2use)^n }

               /    1  (1/2)        1  (1/2)   1\ 
              { C = - 5     , Phi = - 5      + - }
               \    5               2          2/ 

# FORMULA 1

F2 := n -> C * (Phi^n -(1-Phi)^n):

ToProve_1 := F2(n + 1)^2 - F2(n)*F2(n+2) = (-1)^n:

expand~(ToProve_1):
simplify(eval(%, rel2use))
                             n       n
                         (-1)  = (-1) 

Formula 2:

# Preliminary:

G2 := n -> arccot(``(F2)(n));
              G2 := n -> arccot((F2)(n))

ToProve_2 := G2(2*n) - (G2(2*n + 1) + G2(2*n + 2)) = 0;
    arccot((F2)(2 n)) - arccot((F2)(2 n + 1)) - arccot((F2)(2 n + 2)) = 0

# From recurrence formula, replacing n+1 ny 2*n+2 and F by F2:

recur := ``(F2)(2*n+2) = ``(F2)(2*n+1) + ``(F2)(2*n);
           (F2)(2 n + 2) = (F2)(2 n + 1) + (F2)(2 n)

# One gets:

subs(recur, ToProve_2);
 arccot((F2)(2 n)) - arccot((F2)(2 n + 1)) - arccot((F2)(2 n + 1) + (F2)(2 n)) = 0

# Which has symbolic form:
# (U > 0, V > 0, U < V)

SymbolicForm := arccot(U) - arccot(V) - arccot(U+V)
             arccot(U) - arccot(V) - arccot(U + V)


# Proof

G2 := n -> arccot(F2(n)):

recur := F2(2*n+2) = F2(2*n+1) + F2(2*n):

ToProve_2 := simplify(G2(2*n) - (G2(2*n + 1) + G2(2*n + 2)) = 0):

SymbolicForm := -(1/2)*Pi + arctan(U) - arctan(V) - arctan(U+V)
           1                                           
         - - Pi + arctan(U) - arctan(V) - arctan(U + V)
           2                                           

# Thus 

SymbolicForm := -(1/2)*Pi + arctan((U-V)/(1+U*V)) - arctan(U+V):

# and next

SymbolicForm := simplify( -(1/2)*Pi + arctan(((U-V)/(1+U*V) - (U+V)) / (1+(U-V)/(1+U*V)*(U+V))) )
                              /  / 2          \ \
                 1            |V \U  + U V + 2/ |
               - - Pi - arctan|-----------------|
                 2            | 2          2    |
                              \U  + U V - V  + 1/

FocusOn := denom( op(1, indets(SymbolicForm, function)[]) )
                        2          2    
                       U  + U V - V  + 1

FocusOn := eval(FocusOn, {U=F2(2*n), V=F2(2*n+1)}):
FocusOn := eval(FocusOn, rel2use):
FocusOn := simplify(FocusOn) assuming n::integer;
                         (4 n)                     (4 n + 1)
      1  /  1  (1/2)   1\        1 /  1  (1/2)   1\         
    - -- |- - 5      + -|      + - |- - 5      + -|         
      10 \  2          2/        5 \  2          2/         

                            (4 n)       
         1  /  1  (1/2)   1\       (1/2)
       + -- |- - 5      + -|      5     
         10 \  2          2/            

FocusOn := ``(Y^(4*n)) * expand(eval(FocusOn, (-(1/2)*sqrt(5)+1/2)=Y) / Y^(4*n))
               / (4 n)\ /  1    1     1   (1/2)\
               \Y     / |- -- + - Y + -- 5     |
                        \  10   5     10       /

eval(FocusOn, Y=-(1/2)*sqrt(5)+1/2);
                               0
# Then 

'SymbolicForm' = op(1, SymbolicForm) + arctan(something/``(0))
                            1            /something\
           SymbolicForm = - - Pi + arctan|---------|
                            2            \   (0)   /

something := numer( op(1, indets(SymbolicForm, function)[]) )
                          / 2          \
                        V \U  + U V + 2/

# Obviously 'something' is strictly positive, thus 
#
# arctan(something/``(0)) = Pi/2;
#
# and thus

'ToProve_2' = ``(-Pi/2) + ``(Pi/2)

                           /  1   \    / 1   \     
               ToProve_2 = |- - Pi| + |  - Pi|
                           \  2   /    \ 2   /     

Download Fibonacci.mw

PS I believe the second proof could be easier to catch if one uses the trigonometric representation of the golden number.

@jalal 

" all calculators are programmed to divide by n... it's a broad debate"
Not all the calculators are programmed to divide by n (look to Excel^(*) which proposes to divide by n or by n-1).
And it is not a broad debate: the situation where you have to divide by n is extremely clear, as is the situation where you have to divide by n-1.
The only situation when you have to divide by n is when your sample represents the whole population. If it is not the case dividing by n lead to a biased estimator of the variance, while dividing by n-1 lead to an unbiased estimator.

The reason I used n-1 (without asking you which one of these two situations you were in) is a direct consequence of the way you construct data and Weights:

• As you randomly chose the elements of data from the range 4..20 this implicitely means that you draw a sample from a population whose indicuals are 4, 5, ..., 20. So data does not represents the whole population.
• The same reasoning applies to Weights which are randomly chosen.

Had you said "I have population made of W[1] individuals X[1], ..., W[N] individuals X[N]" , the situation would have been different and the division by n perfectly justified.

Excel^(*) : Excerpt from  Variance Excel:

               Excerpt from Population variance Excel

REMARK: Alternative ways to generate data are

N := 10:

data := sort(combinat:-randperm([$4..20])[1..N]);
data := sort(Statistics:-Shuffle([$2..20])[1..N]);


CodeTools:-Usage( sort(Statistics:-Shuffle([$2..20])[1..N]), iterations=10^5):
memory used=4.73KiB, alloc change=0 bytes, cpu time=46.96us, real time=47.03us, gc time=2.25us

CodeTools:-Usage( sort(combinat:-randperm([$4..20])[1..N]), iterations=10^5):
memory used=2.48KiB, alloc change=0 bytes, cpu time=15.79us, real time=15.80us, gc time=1.14us

@jalal 

...that you have weighted data and that you want to draw a boxplot of these data?

if it is so the trick is to build an unweighted sample from the data and their weights and draw the corresponding boxplot.
This is likely a hammer to kill a fly, but the method which consists in coding your own bowplot is far more complex (see Add-on 2 in attached file number two)


If the weights are strictly positive integers (which is the simpler case) use this:
 

Download Boxplot_of_weighted_data.mw



If weights are strictly positive real numbers (a little bit more complex case) use this file
(look more specifically  the Maple issue of Add-on 3).

 

Download Boxplot_of_real_weighted_data.mw

NOTE:
There exist several quantile estimators for weighted data and the one I used here is just the most commonly used and easier to understand.
See for instance Akinshin

To moderators: please stop converting this comment into an answer, thank you


the expression you want to textplot is the result of the optimization of a Lenard-Jones model you use in your previous question?

if it is so the attached file contains a full reply (likely not the same data points the LJ model comes from).
 

Download display_formula_mmcdara.mw

(UPDATED: please do not turn this comment into an answer)


In your procedure the same quantity m is computed each time the loop is executed. So when theloop ends m contains the last result (here 11).
If you want to procedure to return all the values of m, you must record them in some structure.

For instance (tst returns a list)

tst := proc (M) 
  local i, m; 
  m := NULL; 
  for i from 2 to numelems(M)+1 do 
    m := m, M[1]-M[i-1] 
  end do; 
  return [m] 
end proc:

Or (tst returns a vector)

tst := proc (M) 
  local i, m; 
  m := Vector(numelems(M)):
  for i from 2 to numelems(M)+1 do 
    m[i-1] := M[1]-M[i-1] 
  end do; 
  return m
end proc:

Or more simply;

M[1] -~ M[2..-1]

@mz6687 

Here is the solution but it doesn't look to the image you sent before:
 

Download solA_mmcdara.mw


If you really want this "imaginary" form your imagedisplays (I don't find any interest in it), you can obtain it after some dymnastis:

eval(lhs(R), op(1, lhs(R)) = -op(1, lhs(R))*J^2):
newR := op(2, (combine(simplify(%), radical) assuming J > 0)) = Y:
eval( collect(subs(R, valsol[2][1]), Z), Z=J*Y):

s, r := selectremove(has, [op(rhs(%))], J):
convert(add(s) + J^2*(-add(r)), horner, J):
J * value(``(numer(op(1, %))) / denom(op(1, %))):

lambda = subs(J=I, %);
With = ( (rhs=lhs)(newR) )
                        1                        
               lambda = - I (I a[1] + I a[2] + Y)
                        4                        
         /                                             (1/2)\
         |    /    2         2                       2\     |
  With = \Y = \12 c  + 3 a[1]  - 6 a[1] a[2] + 3 a[2] /     /

@Andiguys 

# First thing: construct TRC correctly

TRC_ok := eval(TRC, {alpha = (u -> 0.02*exp(0.4*u)), beta = (u -> 0.2e-1*(1-exp(.5*I2)))})
                        12                   
          4.267000015 10   + 90 Q2 - 21.25 Q1

                             7               
             + 4.000000000 10  I1 exp(0.4 I1)

             + 2000000000 I2 (0.02 - 0.02 exp(0.5 I2))

# Assuming Q2 is nonnegative, TRC_ok is obviously mimimum when Q2=0.

TRC_reduced := eval(TRC_ok, Q2=0);
                12                            7               
  4.267000015 10   - 21.25 Q1 + 4.000000000 10  I1 exp(0.4 I1)

     + 2000000000 I2 (0.02 - 0.02 exp(0.5 I2))

# To be sure of what you are about to do verify what are the indeterminates of TRC_ok:

indets(TRC_ok, name)
                        {I1, I2, Q1, Q2}

# As you commented C1, C2 and C3, the Minimize command cannot work

C1 := 4e8 >= Q1,  Q1 >= .6e7;
C3 := I1 <= I2;
                              8      6      
                    Q1 <= 4 10 , 6 10  <= Q1
                            I1 <= I2

Next:

# Do you really want to write this:
#  C2 := if Q1= 4*(10)^(8) then Q2 =0 else Q2 = (x*d*h)-(delta*Q1*h)
# which meand that Q2 is always equal to (x*d-delta*Q1)*h excepted when Q2=4e8 where Q1=0 ???
#
# As you plane to use a numerical algorithm (Minimize) the probability that Q2=4e8 at some 
# iteration is exactly null... unless the minimum verifies Q1=4e-8 (limit of the first C1 
# condition).

Can we find "by hand" the minimum you are looking for?
TRC_reduced is obviously minimal when Q1 is maximal (Q1=4e8).

# The expression of TRC_reduced shows it is minimal when Q1 is maximal, thus 

evalf(eval(TRC_reduced, Q1=4e8));
               12         7                        9
       4.258 10   + 4.0 10  I1 exp(0.4 I1) + 2.0 10  I2 (0.02 - 0.02 exp(0.5 I2))

As a simple plot shows, the derivative wrt I2 of this expression is negative for 0 <= I2 <= 6.
Thus the minimum of TRC_reduced corresponds to Q1=4e8 and  I2=6.
It remains 

evalf(eval(TRC_reduced, [Q1=4e8, I2=6]));
                        12         7               
                4.254 10   + 4.0 10  I1 exp(0.4 I1)

... which is obviously minimum for I1=0 (a value which verifies I1 <= I2):

evalf(eval(TRC_reduced, [Q1=4e8, I2=6, I1=0]));
                             12
                     4.254 10  

(By the way I think you already asked  this same question and I already answered it. As I don't have time to verify this, please forgive me I'm wrong).

Whatever, here is a way to build an histogram from couples (bins, counts).
I assume your data are a list of equalities bin=count, where bin=a..b is a range and count is an integer.
For other structures, for instance a two-column matrix whose column 1 contains the bins and column 2 the counts, proceed as explained at the end oftheattached file.

Bins_and_counts_histogram.mw

I remain at your disposal for further help.

Remark: in the arrached file STEP 2 shows how to get the bins and counts of raw data while using TallyInto.

If you want to estimate the size and shape parameters of the Weibull random variable the data S could bea sample drawn from, you just have to do this

mle := MaximumLikelihoodEstimate(Weibull(a, b), S);
                 mle := [a = 129.944334168688, b = .9301095062]

If you want to build the corresponding random variable for any other purpose, I advice you to do this

W   := RandomVariable(Weibull(a, b)):
mle := MaximumLikelihoodEstimate(W, S):
W__opt := Specialize(W, mle):

# Example of use

N := numelems(S);
QuantilePlot(Sample(W__opt, N), S):

Point 1: Your differential system contains 8 parameters you have to fix before solving it numerically 

ind := convert(indets(ODE union IC, name) minus {eta}, list);
         [Le, Pr, alpha, N[b], N[t], a[1], a[2], a[3]]

Point 2: When you write

dsolve(eval(ODE union IC, a[i]=a), ...

what index "i" are you refering to?

Point 3: the output of your procedure finda is 

eval(f(eta),g(eta),h(eta),tt(10));

But eval needs two and only two argiments: what does this command of yours is meant to do?

Point 4: if you meant this

eval([f(eta),g(eta),h(eta)],tt(10));

the output of finda is a list with 3 elements, so 

fsolve(finda(a)-1, a=0..1);

is a nonsense because a list cannot be equal to a scalar ("1").

Look here for more explanations Shoot_Blasius-mmcdara.mw

See help(LinearAlgebra[Eigenvectors]): the output of Eigenvectors is a sequence val, vec wher val is a vector containing the eigenvalues and vec the matrix of eigenvectors.
To simplify the output of Eigenvectors do this

simplify([Eigenvectors(L)])

# or

val, vec := Eigenvectors(L):
simplify(vec)

Here is a more synthetic expression of the eigenvalues and eigenvectors:
synthetic.mw

Be aware of undesired effects produced by break, look to this example

restart
for n from 1 to 4 do
  for i from 1 to 2 do 
    if n < 2 then 
      printf("   Within if structure %d  %d\n", n, i)
    else
      break
    end if; 
    printf("Exit if structure %d  %d\n\n", n, i)
  end do;
end do;

   Within if structure 1  1
Exit if structure 1  1

   Within if structure 1  2
Exit if structure 1  2

So just replace

        if Norm(err) <= tol then 
            h := 0.1 * h * (c + 1) * (tol/Norm(err))^(0.2);
        else 
            break
        end if;

by

        if Norm(err) <= tol then 
            h := 0.1 * h * (c + 1) * (tol/Norm(err))^(0.2);
        end if;

Two comments:

Are you sure Digits:=30 is necessary and that  Digits:=10 (default value) would npt ne enough?

Given the large number of points I would prefer one of these two graphics (but it is personal opinion):

numerical_plot_y1 := plot(time_t, numerical_array_y1, style = point, symbol=point, color = blue, legend = ["TFIBF"]):
exact_plot_y1 := plot(time_t, exact_array_y1, style = point, symbol=point, color = red, legend = ["EXACT"]):
display({numerical_plot_y1, exact_plot_y1}, size=[800, 400]);

# or

numerical_plot_y1 := plot(time_t, numerical_array_y1, style = line, color = blue, legend = ["TFIBF"]):
exact_plot_y1 := plot(time_t[1..-2], exact_array_y1[1..-2], style = line, color = red, legend = ["EXACT"]):
display({numerical_plot_y1, exact_plot_y1}, size=[800, 400]);

In Maple 2015 (the version I use) printlevel is an argument of kernelopts, not of interface.
Did you check that?
Note that you can also write printlevel := n directly.

Just a question: Why do you want to change the value of printlevel within your procedure?
If you want to display intermediate results I recommend you to use trace instead:

restart
fac := proc(n::integer)
    if n = 0 then
        1
    else
        n * thisproc(n - 1);  
    end if;
end proc;
proc(n::integer)  ...  end;
fac(5)
                              120
trace(fac):
fac(5)

{--> enter fac, args = 5
{--> enter fac, args = 4
{--> enter fac, args = 3
{--> enter fac, args = 2
{--> enter fac, args = 1
{--> enter fac, args = 0
                               1
<-- exit fac (now in fac) = 1}
                               1
<-- exit fac (now in fac) = 1}
                               2
<-- exit fac (now in fac) = 2}
                               6
<-- exit fac (now in fac) = 6}
                               24
<-- exit fac (now in fac) = 24}
                              120
<-- exit fac (now at top level) = 120}
                              120

Next proceed as described in the attached file:
Ode_New_TWO_phase_mmcdara.mw

Advice...
Odes 2 and 4 both form a subsystem with unknown functions h(Theta) and H(Teta) that you should solve apart.
Note that it can be solved numerically or formally.
Once done, plug this soution into odes 1 and 3 (it's better to have solved the first substystem formally, but you can do the same using the option known [see help(dsolve[numeric]) ] if you solved the first subsystem numerically).

Errors to fix: 

• The first subsystem requires 3 ic/bc and you give 4.
• The first subsystem requires 5 ic/bc and you give 4.

@Prakash J 

Look at this file Difference_mmcdara.mw

@C_R Feel free to delete or reappropriate this comment if you consider that I have interfered inappropriately in your exchanges (I am no longer interested in this race for badges which seems to me to pollute Mapleprimes :-) )

To moderators: this post is deliberately a comment and not an answer, so please don't convert it into an answer. Thanks.