I have constructed a formula for the terms of the Fibonacci sequence for positive integer $n$. It is: $$ \frac{a^{n}}{1+b^2}+\frac{b^n}{1+a^2} $$ where $a=\frac{1+\sqrt{5}}{2}$, and $b=\frac{1-\sqrt{5}}{2}$. I tested few terms out and it does represent the Fib. sequence well, but I want to have a concrete mathematical proof of it as well. Which I planned to do by referring to a fact that the division of two neighboring Fib. terms approaches the Golden ratio, $\frac{1+\sqrt{5}}{2}=a$ as $n$ approaches $\infty$. Therefore I wanted to prove my formula for Fib. terms by showing that
$$ \frac{\frac{a^{n+1}}{1+b^2}+\frac{b^{n+1}}{1+a^2}}{\frac{a^{n}}{1+b^2}+\frac{b^n}{1+a^2}}=a. $$ By taking advantage of the fact that $a=-b^{-1}$ and $b=-a^{-1}$ I arrived to the following $$ \frac{\frac{a^2a^{n+1}+b^{n+1}}{1+a^2}}{\frac{a^2a^n+b^n}{1+a^2}}=\frac{a^2a^{n+1}+b^{n+1}}{a^2a^n+b^n} $$ And from here since $b=-a^{-1}$, $b^{n+1}=\frac{-1}{a^{n+1}}$ and $b^{n}=\frac{-1}{a^{n}}$. When we substitute these values back we get: $$ \frac{a^{n+3}+\frac{-1}{a^{n+1}}}{a^{n+2}+\frac{-1}{a^{n}}}=\frac{\frac{a^{2n+4}-1}{a^{n+1}}}{\frac{a^{2n+3}-1}{a^{n}}}=\frac{a^{2n+4}-1}{a(a^{2n+3}-1)} $$ And at this moment I was not sure how to advance. Any help in proving this fraction is equal to $a$ will be appreciated, thank you.
As @kandb noted in the comments, the correct mathematical statement should say that the limit equals to the golden ratio, not that any adjacent terms divide to the golden ratio. That is, you are saying that you want to prove
$$ \lim_{n \to \infty} \frac{\frac{a^{n + 1}}{1 + b^2} + \frac{b^{n + 1}}{1 + a^2}}{\frac{a^n}{1 + b^2} + \frac{b^n}{1 + a^2}} = a $$
(Bravo for typing it correctly!) Actually, you have come really close already. In your last simplifying line, you wrote
$$ (LHS = ) \lim_{n\to\infty}\frac{a^{n+3}+\frac{-1}{a^{n+1}}}{a^{n+2}+\frac{-1}{a^{n}}}=\lim_{n\to\infty}\frac{\frac{a^{2n+4}-1}{a^{n+1}}}{\frac{a^{2n+3}-1}{a^{n}}}=\lim_{n\to\infty}\frac{a^{2n+4}-1}{a(a^{2n+3}-1)} $$
However, there is an arithmetic mistake:
$$ (LHS = ) \lim_{n\to\infty}\frac{a^{n+3}+\frac{-1}{a^{n+1}}}{a^{n+2}+\frac{-1}{a^{n}}}=\lim_{n\to\infty}\frac{(a^{2n+4}-1)/a^{n+1}}{(a^{2n+{\color{red}2}}-1) / a^{n}}=\lim_{n\to\infty}\frac{a^{2n+4}-1}{a(a^{2n+{\color{red}2}}-1)} $$
Fixing this mistake, the problem becomes easy. It is not hard to show that the constants $- 1$ don't matter, i.e.
$$ \lim_{n\to\infty} \frac{a^{2n + 4} - 1}{a(a^{2n + 2} - 1)} = \lim_{n\to\infty} \frac{a^{2n + 4}}{a(a^{2n + 2})} = a $$
Where the right hand side is exactly what you want to prove! This "removing constants" can be justified by various ways e.g. the squeeze theorem. I can elaborate on this if you need it :)
Just for your interest (and @kandb's), I will show that all second-order linear recurrences with positive coefficients will grow. Specifically, consider a recurrence $f(n) = af(n - 1) + bf(n - 2)$ and some initial values $f(1) = x$ and $f(2) = y$ that won't matter. Notice(!) that we can write this as a matrix equation:
$$ \begin{pmatrix} f(n) \\ f(n - 1) \end{pmatrix} = \begin{pmatrix} a & b \\ 1 & 0 \end{pmatrix} \begin{pmatrix} f(n - 1) \\ f(n - 2) \end{pmatrix} $$
By some simple induction, we see that
$$ \begin{pmatrix} f(n) \\ f(n - 1) \end{pmatrix} = \begin{pmatrix} a & b \\ 1 & 0 \end{pmatrix}^{n - 2} \begin{pmatrix} f(2) \\ f(1) \end{pmatrix} = A^{n - 2}\vec{v} $$
Now this is just a standard linear algebra problem. If you work out all the details, diagonalising the matrix and such, you will see that you get $f(n) = A\lambda_1^n + B\lambda_2^n$ for some unspecified constants $A, B$, where $\lambda_1$ and $\lambda_2$ are roots of the equation $\lambda^2 = a\lambda + b$. Look - this looks exactly like our recurrence definition $f(n) = af(n - 1) + bf(n - 2)$! Solving this equation for $a = b = 1$ also yields $\lambda = \frac{1 \pm \sqrt{5}}{2}$, which appears in Pinet's formula that you wrote. This can be shown via induction, or Cayley-Hamilton, or whatever. Also, this equation $\lambda^2 = a\lambda + b$ is called the characteristics polynomial of the recurrence, which generalises naturally to higher order recurrence.
Finally, to link this back to your original question of finding the limit of the ratio of adjacent terms, we look at
$$ \lim_{n \to \infty} \frac{A\lambda_1^{n + 1} + B\lambda_2^{n + 1}}{A\lambda_1^n + B\lambda_2^n} $$
I forgot about the details exactly, but for "nice" values of $\lambda$'s like $\max(\lambda_1, \lambda_2) > 1$ and $\min(\lambda_1, \lambda_2) > -1$, it is easy to show that the limit is $\max(\lambda_1, \lambda_2)$. Intuitively it makes sense: In the long run, the expression $A\lambda_1^n + B\lambda_2^n$ is dominated by whichever $\lambda$ is larger, modulo some negative annoyance.
Here's some extra questions to think about. You can also construct a sequence with the properties below and look at the values numerically :)