Let $f_n: \Bbb N_0 \rightarrow \Bbb N_0$ be a function for every $n\in \Bbb N_0$. Let $(f_n)_{n=1}^\infty$ be a sequence given by the following recurrence-relation:
$$ f_1(k) = 1 $$ $$ f_2(k) = k $$ $$ f_n(k + 1) = f_n(k) + f_{n-1}(k) + f_{n-2}(k) + f_{n-1}(k)f_{n-2}(k)$$
Futhermore it is assumed that $f_n(0) = 0$ for every $n > 2$.
I am able to work out what each function would be by hand. The next two functions of the sequence are:
$$ f_3(k) = k^2 $$ $$ f_4(k) = \sum_{i=0}^{k-1} i + i^2 + i^3 = \frac{1}{12} k(k-1)(3k^2 + 4k + 4)$$
I am trying to find a general formula for what $f_n$ would look like for arbitrary $n$, but i don't know how to proceed. Does such a formula exist and how would i go about finding it?
TLDR summary so far: we know the integer roots of these $f_n$, and we know that the sum of their coefficients is always zero ($n\ge4$) (and the constant coefficient is zero for all $f_n$, $n\neq1$). We also have an ugly formula for the degree of the polynomial for each $f_n$, and a way of computing the highest order coefficients of $f_n$.
Ongoing Answer:
Each of the $f_n$ is a polynomial of degree $d$: $a_0+a_1k+a_2k^2+\cdots+a_dk^d$ and it can be seen that $f_n(k+1)-f_n(k)$ gives a polynomial of degree $d-1$. The multiplication of two polynomials of degree $a$ and $b$ gives a polynomial of degree $a+b$. Their sum gives a polynomial of degree $\max(a,b)$.
Rearrange your recurrence relationship to: $f_n(k+1)-f_n(k)=f_{n-1}(k)+f_{n-2}(k)+f_{n-1}(k)\cdot f_{n-2}(k)$. Take the degree of both sides:
$$\deg [f_n(k+1)-f_n(k)]=\deg[f_{n-1}]+\deg[f_{n-2}]\implies d-1=\deg[f_{n-1}]+\deg[f_{n-2}]$$
Reconstructing from the base cases of $\deg[f_1]=0,\deg[f_2]=1$, we get something like Fibonacci:
$$\deg[f_1]=0\\\deg[f_2]=1\\\deg[f_3]=1+\deg[f_1]+\deg[f_2]=1+0+1=2\\\deg[f_4]=1+2+1=4\\\deg[f_5]=1+2+4=7\\\cdots$$
From hereon out I will call "$\deg[f_n]$" $D_n$.
Sometimes with linear recurrence relationships, we can set it up in a matrix form.
$$\begin{bmatrix}D_n\\D_{n-1}\\D_2=1\end{bmatrix}=\begin{bmatrix}1&1&1\\1&0&0\\0&0&1\end{bmatrix}\begin{bmatrix}D_{n-1}\\D_{n-2}\\D_2=1\end{bmatrix}\implies\begin{bmatrix}D_n\\D_{n-1}\\D_2=1\end{bmatrix}=\begin{bmatrix}1&1&1\\1&0&0\\0&0&1\end{bmatrix}^{n-2}\begin{bmatrix}D_2=1\\D_1=0\\D_2=1\end{bmatrix}$$
You compute the $n$th power of a matrix by diagonalising (when possible, which it is here). I will let $\phi=\frac{1}{2}(1+\sqrt{5})$ and $\varphi=\frac{1}{2}(1-\sqrt{5})$ - a true similarity with the Fibonacci series. Anyway, the diagonalisation $M=SJS^{-1}$ looks like this:
$$\begin{bmatrix}1&1&1\\1&0&0\\0&0&1\end{bmatrix}^{n-2}=\begin{bmatrix}-1&\varphi&\phi\\-1&1&1\\1&0&0\end{bmatrix}\begin{bmatrix}1&0&0\\0&\varphi^{n-2}&0\\0&0&\phi^{n-2}\end{bmatrix}\begin{bmatrix}0&0&1\\-\frac{1}{\sqrt{5}}&\frac{\phi+2}{5}&\frac{\varphi+2}{5}\\\frac{1}{\sqrt{5}}&\frac{\varphi+2}{5}&\frac{\phi+2}{5}\end{bmatrix}$$
And we only care about $D_n$, so only the top row of the matrix will be considered. That row comes out as:
$$\left(\frac{\phi^{n-1}-\varphi^{n-1}}{\sqrt{5}},\frac{\varphi^{n-1}(\phi+2)+\phi^{n-1}(\varphi+2)}{5},\frac{\varphi^n+\phi^n+2(\varphi^{n-1}+\phi^{n-1})-5}{5}\right)$$
Ugly, right? If we keep pushing on, we will get a formula for the degree of the $n$th polynomial. Multiply this top row by the initial values vector $(1,0,1)^T$, and we get:
$$\begin{align}D_n&=\frac{\phi^n+\sqrt{5}\cdot\phi^{n-1}+\varphi^n-\sqrt{5}\cdot\varphi^{n-1}+2\phi^{n-1}+2\varphi^{n-1}}{5}-1\\\deg[f_n(k)]&=\frac{\phi^n+\varphi^n+\phi^{n-1}(2+\sqrt{5})+\varphi^{n-1}(2-\sqrt{5})-5}{5}\end{align}$$
Ta-Da! And I have tested this in python - it does work!
The next bit is an induction to figure out the integer roots of these polynomials.
You said that by definition that for all $f_n,n\gt2$, $f_n(0)=0$. We can actually do better than that. Look at the recurrence relation: if $k$ is a root of $f_{n-1}$, $f_{n-2}$, and a root of $f_n$, then $k+1$ is also a root of $f_n$.
$$f_4(1)=f_4(0)+f_3(0)+f_2(0)+f_3(0)\cdot f_2(0)=0\\f_5(1)=\cdots=0\\f_6(1)=0\\\vdots\\f_6(2)=f_6(1)+f_5(1)+f_4(1)+f_5(1)\cdot f_4(1)=0\\\vdots\\f_8(3)=0$$
And in fact this statement holds true: for all $n\ge2(k+1)$, $f_n(k)=0$. This tell us that for some $f_n$, we know the first $\lfloor\frac{n}{2}\rfloor$ roots. By the fundamental theorem of algebra, a complex polynomial of degree $n$ has $n$ roots. These are natural number polynomials, but we can still be absolutely sure that for some $f_n$ of degree $D_n$, if we can find $r\lt D_n$ distinct roots, then there must be more to the polynomial. That made little sense, so as an example, take $f_4$. $4=2\cdot2=2(1+1)$ therefore as I just showed there must be roots at $0$ and also at $1$. Since $f_4(1)=0,f_4(0)=0$, by the factor theorem (and also by what I just said) $f_4(k)=k(k-1)\cdot g(k)$, where $g(k)$ is a polynomial of degree $D_4-2=2$ that accounts for the rest of the ($2$ non-natural) roots. Remember all non-constant polynomials can look like this: $f(x)=a_d(x-r_1)(x-r_2)(x-r_3)\cdots$ where the $r_i$ are roots and $d$ is the degree. Another way of characterising polynomials is this: $f(x)=a_0+a_1x+a_2x^2+\cdots$ Our recurrence relation tells us $f_n(k+1)-f(k)$ in terms of the lower-order polynomials. $f_n(k+1)-f(k)$ has a lower degree, as I said before, but I didn't prove that. I'll show it now, alongside an analysis of the expression.
$$\begin{align}f_n(k+1)-f_n(k)&=a_0+a_1(k+1)+a_2(k+1)^2+\cdots-a_0-a_1k-a_2k^2-\cdots\\&=\sum_{i=1}^{D_n}a_i\sum_{j=0}^{i-1}\binom{i}{j}k^j=a_1+a_2(2k+1)+\cdots+a_{D_n}(k^{D_n-1}+\cdots)\\&=\sum_{i=0}^{D_n-1}k^i\sum_{j=i+1}^{D_n}\binom{j}{i}a_j=(a_1+a_2+\cdots)+k(2a_2+3a_3+\cdots)+k^2(3a_3+6a_4+\cdots)+\cdots\\&=f_{n-1}(k)+f_{n-2}(k)+f_{n-1}(k)\cdot f_{n-2}(k)\\&=(b_1k+b_2k^2+\cdots)+(c_1k+c_2k^2+\cdots)+(b_1k+b_2k^2+\cdots)(c_1k+c_2k^2+\cdots)\end{align}$$
The second line shows you why $\deg[f_n(k+1)-f_n(k)]=D_n-1$. The second line also amounts to the sum of all the coefficients (except $a_0$, but we know $a_0$ equals zero anyway since $0$ is a root to all the polynomials $n\ge2$), when $k=0$. But we know that $f(1)$ is often equal to zero! In fact, $f_n(1)=0,\,\forall n\ge4$. So, for all $n\ge4$ and with $k=0$, $f_n(k+1)-f_n(k)=\sum_{i=1}^{D_n}a_i=0$. So, the sum of all the coefficients must equal zero, and $a_0$ is zero, for $n\ge4$. $a_{D_n}$ is also non-zero, since the polynomials have degree $D_n$. We can write:
$$\frac{\sum_{i=0}^{D_n}a_i}{a_{D_n}}=0$$
But Vieta's formulas are interesting here.
$$\begin{align}\frac{\sum_{i=0}^{D_n}a_i}{a_{D_n}}&=0\\&=\frac{a_{D_n}}{a_{D_n}}+\frac{a_{D_n-1}}{a_{D_n}}+\frac{a_{D_n-2}}{a_{D_n}}+\cdots\\&=1-(r_1+r_2+r_3+\cdots)+(r_1r_2+r_1r_3+\cdots+r_2r_3+r_2r_4+\cdots)-\cdots\\&=1-r_1(1-(r_2+r_3+\cdots)+(r_2r_3+r_3r_4+\cdots))\\&-r_2(1-(r_3+r_4+\cdots)+(r_3r_4+r_3r_5+\cdots))-r_3(\cdots)-\cdots\end{align}$$
And I would like to believe that there is some recursion going on here, that I, you, or some reader might be able to build on.
I'm not sure what to make of the recursive Vieta's formula pattern there; but the formulas do tell us something concrete, since $r_1$ is always zero, $n\ge2$.
$$(-1)^{D_n-1}\frac{a_1}{a_{D_n}}=\prod_{i=2}^{D_n}r_i=p!\prod_{i=p+2}^{D_n}r_i$$
I will show you how I calculated the real formula for $f_4$.
By some of the rambling above, we know that for $f_4=a_1k+a_2k^2+a_3k^3+a_4k^4$, $\Delta f_4=(a_1+a_2+a_3+a_4)+k(2a_2+3a_3+4a_4)+k^2(3a_3+6a_4)+k^3(4a_4)=k+k^2+k^3$, by comparing the expansion with the recurrence construction from $f_2,f_3$ which we already know. Comparing coefficients cleanly and unambiguously gives $a_4=\frac{1}{4},a_3=-\frac{1}{6},a_2=\frac{1}{4},a_1=-\frac{1}{3}$ and we know that always $a_0=0$. Therefore:
$$f_4(k)=-\frac{1}{3}k+\frac{1}{4}k^2-\frac{1}{6}k^3+\frac{1}{4}k^4$$
Which satisfies all the requirements stated: it has roots at $0$ and $1$, its coefficients sum to zero, and its degree is the degree given by my formula for $D_n$. The disadvantage of this method is that you need to already know $f_3,f_2$ to solve for $f_4$. Maybe there's a better way, maybe there isn't. I'll keep looking for a while yet.
The polynomial remainder theorem tells us that for $2p+4\ge n\ge2p+2$, $f_n(k)=k(k-1)\cdots(k-p+1)(k-p)Q(k)$. From now on, $q_i$ denotes a coefficient of order $i$ in the quotient polynomial $Q$, and a subscript $a_D$ refers to the highest order coefficient of any polynomial that $a$ refers to - not to be confused with the $D_n$ as a formula for the degree of $f_n$. The superscript $^n$ refers to the $n$th polynomial, for example, $q^n_i$ refers to the $i$th coefficient of the quotient polynomial for $f_n$. So, we know $\Delta f_n=f_{n-1}+f_{n-2}+f_{n-1}f_{n-2}$. We also know that $f_{n-1},f_{n-2}=k(k-1)\cdots(k-p+1)Q^{n-1},Q^{n-2}$ respectively (the bracketed terms are shared). The highest order coefficient of $\Delta f_n$ we know to be $a_{D_n}D_n$, where $a$ refers to the coefficients of $f_n$, and is also, by analysis of the product $f_{n-1}f_{n-2}$, $q^{n-1}_D\cdot q^{n-2}_D$. So $a_{D_n}D_n=q^{n-1}_D\cdot q^{n-2}_D$. We also know that $a_{D_n}=q^n_D$, by a similar analysis. Therefore $q_D=(q^{n-1}_D\cdot q^{n-2}_D)/(D_n)$ - another recurrence relation! What are the initial values? Well, $q^2_D=1,q^3_D=1$, and so rather pleasingly the numerator is all ones, making the highest order coefficient of the quotient polynomial of $f_n$, and thus the highest order coefficient of $f_n$ itself, the following:
$$a_{D_n}=q^n_D=\prod_{i=4}^n(D_i)^{-1}=5^{n-4}\cdot\left(\prod_{i=4}^n\left(\phi^i+\varphi^i+\phi^{i-1}(2+\sqrt{5})+\varphi^{i-1}(2-\sqrt{5})-5\right)\right)^{-1}$$
A similar coefficient analysis that I will not burden the reader with reading, as it is most ugly, yields this equation, which by virtue of being able to solve for $q^n_D$, is solvable to find $q^n_{D-1}$, for all $n\ge5$:
$$q^{n-2}_D\left(\left(\frac{p-p^2}{2}\right)\cdot q^{n-1}_D+q^{n-1}_{D-1}\right)+q^{n-1}_D\left(\left(\frac{p-p^2}{2}\right)\cdot q^{n-2}_D+q^{n-2}_{D-1}\right)=(D_n-1)\left(q^n_{D-1}-\frac{p^2+p}{2}\cdot q^n_D\right)+\frac{D^2_n-D_n}{2}\cdot q^n_D$$
Knowing $q^n_{D_n-1}$ gives us $a_{D_n-1}$, and consequently the sum of all roots, by the formulae:
$$a_{D_n-1}=(D_n-1)\left(q^n_{D-1}-q^n_D\cdot\frac{p^2+p}{2}\right)\\\sum_{i=1}^{D_n}r_i=-\frac{a_{D_n-1}}{a_{D_n}}$$
The subscript "$_{D-1}$", like the subscript "$_D$", I use to refer to the second highest order coefficient of any polynomial. I am using this odd notation because all the polynomials here have different degrees, a nightmare to keep track of!
Once again, the reader is invited to build on my work. For whatever reason, I really like this puzzle and the fact that some progress can be made.