How does one find a formula for the recurrence relation $a_{1}=1,a_{2}=3, a_{n+2}=a_{n+1}+a_{n}?$

There is a standard method which is sketched below. Forget about the initial conditions $a_1=1$, $a_2=3$ for now. We look for solutions of the recurrence relation of the shape $a_n=x^n$, where $x$ is a number to be determined.

Substituting in the recurrence, we get $x^{n+2}=x^{n+1} +x^n$. Forgetting about the possibility $x=0$, which mostly can't happen, this reduces to $x^2-x-1=0$. Solve. We get $$x=\frac{1\pm\sqrt{5}}{2}.$$ Call the two roots $\alpha$ and $\beta$. (The number $\frac{1+\sqrt{5}}{2}$ is a famous number, often called the Golden Number. It even has two standard symbols attached to it, $\varphi$ and $\tau$.)

Now look for numbers $A$ and $B$ such that the expression $A\alpha^n+B\beta^n$ satisfies your initial conditions. (If you have studied Differential Equations, the procedure will be structurally familiar, for good reason.)

The suitable $A$ and $B$ are not hard to find. We get $A=B=1$. To check that $a_n=\alpha^n +\beta^n$ is indeed the solution, note that $A\alpha^n+B\beta^n$ satisfies our recurrence for any $A$ and $B$. Choosing $A$ and $B$ so that the initial conditions are satisfied forces the formula to be correct for all $n$.

It would have been a little more pleasant to find $A$ and $B$ if we had been given the usual initial conditions $a_0=2$, $a_1=1$. Whether we start at $a_0$ or $a_1$, the sequence is called the Lucas sequence.

Remark: Variants of this idea work for all linear homogeneous recurrences with constant coefficients. Look in particular at the recurrence $a_{n+1}=3a_n-2a_{n-1}$ that you solved successfully by thinking.

Use the same procedure. We arrive at the equation $x^2-3x+2=0$, which has roots $\alpha=2$, $\beta=1$. Now we try to find $A$ and $B$ such that $A\alpha^n+B\beta^n$ satisfies our initial conditions. So we want $A(2)+B(1)=3$, $A(2^2)+B(1^2)=5$. We get $A=1$, $B=1$. This yields the formula that you already know.

There are other general methods of doing the same thing, such as generating functions.

Besides Andre Nicolas approach you could also try the following: Generating functions, define (just for a little cleaner solution) $b_0=a_1$, $b_1=a_2$, and in general $b_m=a_{m+1}$.

Let $F(t)=\sum_{n=1}^\infty b_nt^n$. Then $$b_{n+2}=3b_{n+1}-2b_n$$$$b_{n+2}t^{n+2}=3b_{n+1}t^{n+2}-2b_nt^{n+2}$$Taking summation from $n=0$ to infinity: $$\sum_{i=2}^\infty b_it^i=3t\sum_{i=1}^\infty b_it^i-2t^2\sum_{i=0}^\infty b_it^i$$$$F(t)-b_1t-b_0=3t(F(t)-b_0)-2t^2F(t)$$Using the inital conditions:$$F(t)-5t-3=3t(F(t)-3)-2t^2F(t)$$$$F(t)(1-3t+2t^2)=3-4t$$Thus, $$F(t)=\frac{3-4t}{1-3t+2t^2}$$ Break it $$F(t)=\frac{2}{1-2t}+\frac{1}{1-t}$$and that $\sum x^k=\frac{1}{1-x}$ where the sum goes form $k=0$ to infinity we get:$$F(t)=2\sum (2t)^k+\sum t^k$$Thus $b_n=2(2^n)+1$

And moving it back to your initial sequence of $a_n$ you get indeed: $$\boxed{a_n=2^n+1}$$

Another standard method is the matrix method, which is the same as the Fibonacci sequence but with different starting vector. The same solution applies.