Note that $\left\{\binom{x-1}{k}\right\}_{k\ge 0}$ forms a basis for $\mathbb{R}[x]$. Then write $p=a_0\binom{x-1}{0}+...+a_d\binom{x-1}{d}$, where $a_d$ is nonzero ($d$ is degree of $p$). Note then that $p(1)+p(2)+...+p(n)=a_0\binom{n}{1}+...+a_d\binom{n}{d+1}$ by the hockey-stick identity. Call the polynomial $a_0\binom{x}{1}+...+a_d\binom{x}{d+1}$ as $r(x)$. Then as polynomials, $r(x)=p(x)q(x)$. $r$ is a degree $d+1$ polynomial in $x$, so $q$ must have degree $1$. Furthermore this reveals that the leading coefficient of $p(1)+...+p(n)$ when written as a polynomial of $n$ in the standard basis is $\frac{a_d}{d+1}$, so $q(x)$ is of the form $\frac{x}{d+1}+\lambda$ for some $\lambda\in \mathbb{R}$. Note that the map $p\mapsto r$ and the map $p\mapsto \frac{x}{d+1}p$ are both linear, and thus the map $p\mapsto r-\frac{x}{d+1}p$ is a linear map from polynomials of degree at most $d$ to polynomials of degree at most $d$, and all $p$ are eigenvectors of this transformation with eigenvalue $\lambda$. Notice that $\binom{x-1}{d},\binom{x}{d},...,\binom{x+d-1}{d}$ are eigenvectors of this map corresponding to $\lambda=0,\frac{1}{d+1},...,\frac{d}{d+1}$. Since they correspond to $d+1$ different eigenvalues (and the space is $d+1$ dimensional), all eigenspaces are $1$-dimensional, and thus all eigenvectors (and thus valid $p$) are multiples of these. It is not hard to see that all multiples of these work as $p$ by the hockeystick formula, so the solutions for $p$ are precisely $p=c\binom{x+k}{d}$ for $d$ a nonnegative integer and $k$ an integer satisfying $-1\le k\le d-1$.

Let the degree of $p$ be $d$; then by the Integral Test, $q$ is linear with leading coefficient $\frac{1}{d+1}$. For all $n\in \mathbb{N}$, we also have \[p(n)q(n)=p(n+1)(q(n+1)-1)\] so it must be true for all complex $x$ as well; and thus $p(x), p(x+1)$ differ in at most one root, counting multiplicity. It follows that $p$ has roots of the form $r,r+1, \ldots , r+d-1$. Now taking $n=1$ gives $p(1)=0$ or $q(1)=1$. The former implies that $r\in \{-d+2, -d+3, \ldots , 1\}$, and we can find the corresponding $q$ for which the identity always holds. The latter implies that $q(x)=\frac{1}{d+1}x+\frac{d}{d+1}$, and you get the solution for $r=-d+1$.

We can also deduce that q is linear by https://en.wikipedia.org/wiki/Faulhaber%27s_formula and finish like the above.

Here is a solution which is completely elementary: Amir Hossein wrote: Find all polynomials $p(x)$ with real coefficients that have the following property: there exists a polynomial $q(x)$ with real coefficients such that $$p(1) + p(2) + p(3) +\dots + p(n) = p(n)q(n)$$for all positive integers $n$. We begin by observing that if $P(x)$ is such a polynomial, then $\frac{P(x)}{c}$ is also such a polynomial. So WLOG we can assume that the leading co-efficient of $P$ is $1$. This implies that the LHS is positive for large enough $n$ which further implies that the leading co-efficient of $Q$ is positive. Now we can rewrite the given equation as $$P(n-1)Q(n-1)=P(n)(Q(n)-1)$$Observe that if $\deg(q)>2$, for large enough $n$ we would have $P(n)>P(n-1)$ and $Q(n)>Q(n-1)+1$. This implies $Q(x)=ax+b$, where $1>a\ge 0$. Now if $a=0$, it's easy to see that $P\equiv 0$. So let's assume that $Q=ax+b$ with $1>a>0$. Hence we have $$P(n-1)(a(n-1)+b)=P(n)(an+b-1)$$Since the LHS and RHS are also polynomials in $n$, we can conclude that $$ P(x-1)(a(x-1)+b)=P(x)(ax+b-1)\qquad(\bullet)$$ We now claim that all the roots of $P$ will be real numbers with no multiplicity. Suppose that $z$ is a non-real root of $P$, since $P$ has real coefficients, $\bar{z}$ is also a root of $P$. Plugging in $x=z,\bar z$ we see that $z-1,\bar z -1$ are also roots of the LHS. But since roots of a linear equation are always reals numbers, we must have $z-1,\bar z -1$ are also roots of $P$. Thus $P$ has infinitely many zeroes, hence $P\equiv 0$. Thus we know that all the roots of $P$ are indeed real and they form an arithmetic progression with common difference = 1. Suppose the roots of $P$ are $\alpha, \alpha+1,\ldots, \alpha+n$. So we must have $\alpha-1=\tfrac{1-b}{a}$ and $\alpha+n=\tfrac{a-b}{a}$. Suppose that one of the roots $\alpha+i$ appear with multiplicity. Hence $\alpha+i-1$ should also appear with multiplicity. Continuing we have $\alpha+i-1$ should also appear with multiplicity, but it's the root of a linear equation, so it can't possibly appear with multiplicity. Thus we must have $$P(x)=(x-\alpha)(x-\alpha-1)\ldots(x-\alpha-n)$$Now we have $$P(1)+P(2)+P(3)+\ldots+P(k)=\frac{1}{n+1}[(k-\alpha+1)(k-\alpha)\ldots(k-\alpha-n+1)(k-\alpha-n)-(1-\alpha)(-\alpha)(-1-\alpha)\ldots(-n-\alpha)]$$$$\implies \frac{1}{n+1}[(k-\alpha+1)(k-\alpha)\ldots(k-\alpha-n+1)(k-\alpha-n)-(1-\alpha)(-\alpha)(-1-\alpha)\ldots(-n-\alpha)]=(k-\alpha)(k-\alpha-1)\ldots(k-\alpha-n)Q(k)$$Since $Q$ is a linear polynomial, it's easy to see that the possible values of $\alpha$ are $1,0,-1,\ldots -n$. Hence all the solutions are $P\equiv c$ and $P(x)=c(x-\alpha)(x-\alpha-1)\ldots(x-\alpha-n)$ where $\alpha\in\{1,0,-1,\ldots, -n\}$. It's easy to see that they indeed satisfy the given equation.

First we consider the case that $P(x)$ is constant, i.e. $P(x)\equiv c$ hence we have for every positive integer $n$: $nc=cQ(n)\longrightarrow\forall n\in\mathbb{N} : c(Q(n)-n)=0$ Which leads us to solutions of $Q(x)\equiv x , P(x)\equiv c$ and $P(x)\equiv 0 , Q(x)\in\mathbb{R}[X]$ Thus we can assume that $P$ is not constant, i.e. The leading coefficient of $P$ is nonzero. Write $$P(x)=a_sx^s+a_{s-1}x^{s-1}+\dots$$and$$Q(x)=b_tx^t+b_{t-1}x^{t-1}+\dots$$Note that the given equation can also be rewritten into $$\forall n\geq2: P(n-1)Q(n-1)=P(n)(Q(n)-1)$$Since this holds for infinitely many n, it holds for all Real numbers, i.e. we have a Polynomial equality, thus if we have both $\deg(P)>1$ and $\deg(Q)>1$ by checking The coefficient of $x^{t+s-1}$ we get $$[x^{t+s-1}]=a_sb_{t-1}+a_{s-1}b_t-sa_sb_t-ta_sb_t=a_sb_{t-1}+a_{s-1}b_t$$and we have $(s+t)a_sb_t=0$ but we had assumed $a_s,b_t\ne0 \text{and} n+s>2$, a contradiction. Thus we must have $\deg(P)\leq1$ or $\deg(Q)\leq1$. Now assume that we have $\deg(P)\leq1$, since we know P is not constant, we must have $\deg(p)=1$. Now if $\deg(Q)>1$ checking the coefficient of $x^t$ we get $$[x^t]=a_1b_{t-1}+a_0b_t-ta_1b_t=a_1b_{t-1}+a_0b_t$$once again leading to a contradiction. We conclude that we must have $\deg(Q)\leq1$. Now we'll prove that of each degree $s$ for polynomial $P$ we have a unique form of solution for $P$ and a first degree $Q$. First assume that we have $\deg(Q)=0$ for some $s$, i.e. $Q(x)\equiv T$ then we have $$P(x)T=P(x)(T-1)$$and checking the leading coefficient leads to $P$ being constant, a contradiction, hence we know that if for $P$ of degree $s$ we have a solution, then $Q$ will be of first degree. Let $Q(x)=bx+c$ then $Q(x)-1=bx+c-1$ and $Q(x-1)=b(x-1)+c=bx+(c-b)$ and once again checking the coefficient of $x^s$ we have $$[x^s]=a_s(c-b)+a_{s-1}b-bsa_s=a_s(c-1)+a_{s-1}b$$Thus we get $a_sb(s+1)=a_s$ or equivalently $b=\frac{1}{s+1}$ (since $a_s$ is nonzero by assumption). Let us write $P$ as follows $$P(x)=\alpha(x-r_1)(x-r_2)\dots(x-r_s) , r_1\leq r_2\leq\dots\leq r_s$$Plugging into the equality $P(x)(Q(x)-1)=P(x-1)Q(x-1)$ the values $x-1=r_s$ and $x=r_1$ we get $$Q(r_s+1)=1 , Q(r_1-1)=0 \longrightarrow Q(r_s+1)-Q(r_1-1)=\frac{r_s-r_1+2}{s+1}=1\longrightarrow r_s-r_1=s-1$$also note that since $Q$ is linear, The Equations $Q(x)=1$ and $Q(x)=0$ have exactly one solution, and since we know these solutions are $r_s+1$ and $r_1-1$ respectively, then none of these equations hold for $r\in(r_1-1,r_s+1)$. Specifically we get that $P(r_1)=P(r_1+1)=\dots=P(r_s)=0$ and thus we have that $r_1=r , r_i=r+i-1$ Thus $P(x)=\alpha(x-r)(x-r-1)\dots(x-r-s+1)$ Note that we have $$P(1)=P(1)Q(1)\longrightarrow P(1)=0 \quad\text{or}\quad Q(1)=1$$If we have $P(1)=0$, we get $r+k=1$ with $0\leq k\leq s-1$ hence we have $r\in\{1,0,\dots,-s+2\}$ If we have $Q(1)=1$ we get $r_s+1=1$ thus $r+s-1=0$ leading to $r=-s+1$ Thus we have $$P(x)=\alpha(x-r)(x-r-1)\dots(x-r-s+1) , \alpha\in\mathbb{R} , r\in\{1,0,\dots,-s+1\}$$One can easily check that these do indeed satisfy the conditions and since we had $Q(r-1)=0$ and $b=\frac{1}{s+1}$, $c$ is found and $Q$ is uniquely defined as $Q(x)=\frac{x+1-r}{s+1}$ We are done.

Idea:for $P(x)$ exist $R(x)$ such that $$P(x)=R(x+1)-R(x)$$Bu Lagrange interpoletion.