Quantising action to find evolution of wavefunction I am given the following action functional $$S=\int dt \left[p_1\dot q_1+p_2\dot q_2-\frac{p_1p_2}{m}-\sum_{n=0}^\infty\frac{\partial^nV(q_1)}{\partial q_1^n}q_2^n\right]$$
where $p_i$ is the momentum conjugate to the coordinate $q_i$. I am then asked to quantise this action as a two dimensional system and show that $\psi:=\psi(q_1,p_2,t)=\left\langle q_1,p_2\middle|\psi(t)\right\rangle$ follows time evolution $$\frac{\partial\psi}{\partial t}=-\frac{p_2}{m}\frac{\partial\psi}{\partial q_1}+\sum_{n=0}^\infty(-\hbar^2)^n\frac{\partial^{2n+1}V(q_1)}{\partial q_1^{2n+1}}\frac{\partial^{2n+1}\psi}{\partial p_2^{2n+1}}$$
My attempt:
I first find the classical Hamiltonian $$H=p_1\dot q_1+p_2\dot q_2-L=\frac{p_1p_2}{m}+\sum_{n=0}^\infty\frac{\partial^nV(q_1)}{\partial q_1^n}q_2^n$$
I was thinking of then using $$\frac{\partial\psi}{\partial t}=-\left\{\left\{\psi,H\right\}\right\}$$ where $\{\{\ \}\}$ is the Moyal bracket, but it doesn't really get me anywhere.

UPDATE My bad. In this question I gave the wrong expression for the action, which in fact was $$S=\int dt \left[p_1\dot q_1+p_2\dot q_2-\frac{p_1p_2}{m}-\sum_{n=0}^\infty\frac{1}{(2n+1)!2^n}\frac{\partial^{2n+1}V(q_1)}{\partial q_1^{2n+1}}q_2^{2n+1}\right]$$
I appreciate Cosmas for taking the time to answer the question. However, it turns out that the answer I was looking for follows quite simply from the conventional quantisation of the action above.
 A: It appears you have imposed an extra condition that the wavefunction be real. The "potential term" appears complex, so it should connect the real and imaginary parts of the wavefunction. (See below).
To quantize your classical 
$$H=\frac{p_1p_2}{m}+\sum_{n=0}^\infty\frac{\partial^nV(q_1)}{\partial q_1^n}q_2^n$$
into a $q_1,p_2$ basis, you readily get
$$\frac{\partial\psi}{\partial t}=-\frac{p_2}{m}\frac{\partial\psi}{\partial q_1}+\frac{1}{i\hbar} \sum_{n=0}^\infty(i\hbar)^n\frac{\partial^{n}V(q_1)}{\partial q_1^{n}}\frac{\partial^{n}\psi}{\partial p_2^{n}}.$$
Recall that, in the momentum basis, the coordinate is represented by minus the gradient you are used to from the coordinate basis ( $\hat p=\int dx |x\rangle (-i\hbar \partial_x )\langle x| $ ),
$$
\hat x= \int  dx ~~ |x\rangle x \langle x| = \int  dx dp~ |x\rangle x \langle x|p\rangle \langle p| =  \int  dx dp~ |x\rangle (-i\hbar \partial_p \langle x|p\rangle )  \langle p| \\  
 =  \int  dx dp~ |x\rangle   \langle x|p\rangle (i\hbar \partial_p ) \langle p| =  \int    dp~  |p\rangle i\hbar \partial_p  \langle p| ,
$$
where you observed the integration by parts in p at the break of the lines.
Now, the "potential" term is complex in general, and so connects the real and imaginary parts of your wavefunction. If you are seeking a real wavefunction, then the imaginary part of the above TDSE must vanish separately, and the real part reduces to the projection of the above to just the odd terms, 
$$\frac{\partial\psi}{\partial t}=-\frac{p_2}{m}\frac{\partial\psi}{\partial q_1}+  \sum_{k=0}^\infty(-\hbar^2)^k\frac{\partial^{2k+1}V(q_1)}{\partial q_1^{2k+1}}\frac{\partial^{2k+1}\psi}{\partial p_2^{2k+1}} ~.$$
Since the operators $\partial_{q_1}$ and $\partial_{p_2}$ commute, you could further compact the infinite series to a more  formal fractional expression,
$$
V(q_1) ~\frac{\overleftarrow{\partial }_{q_1}
\overrightarrow{\partial }_{p_2} } {1+\hbar^2 \overleftarrow{\partial }_{q_1}^2 ~
\overrightarrow{\partial }_{p_2}^2} ~\psi (q_1,p_2)   ~~~.
$$

Update As per the question's update, the last two formulas may represent an unrestricted complex TDSE, without reality restrictions. 
