Cauchy Reimann equations-help with formal proof..

I'm working on a problem which involves two harmonic functions $\displaystyle p$ and $\displaystyle \zeta$ defined in a square.

($\displaystyle p_{xx} + p_{yy} =0$ , $\displaystyle \zeta_{xx} + \zeta_{yy} = 0$ everywhere.)

The boundary conditions are $\displaystyle p,\zeta=0$ along $\displaystyle x=0,1$ and $\displaystyle p,\zeta=sin^2(\pi x)$ along $\displaystyle y=1$.

Along $\displaystyle y=0$ the functions satisfy the Cauchy Riemann equations
$\displaystyle p_x=-\zeta_y$
$\displaystyle p_y= \zeta_x$

I know (via some dimensional analysis from the Navier Stokes equations) that because of the condition at $\displaystyle y=0$ these functions must satisfy the Cauchy Riemann equations everywhere, and can show this with an asymtotic expansion from $\displaystyle y=0$:

$\displaystyle p_x = p_x|_{y=0} +y p_{xy}|_{y=0} + \frac{y^2}{2!} p_{xyy}|_{y=0} + ...$
$\displaystyle \zeta_y = \zeta_y|_{y=0} +y \zeta_{yy}|_{y=0} + ... \frac{y^2}{2!} \zeta_{yyy}|_{y=0} +...$

But, from the C-R at $\displaystyle y=0$,

$\displaystyle p_{xy}|_{y=0} = \zeta_{xx}|_{y=0}$
$\displaystyle p_{xyy}|_{y=0}= \zeta_{xxy}|_{y=0}$

and so on..

So that
p$\displaystyle _x + \zeta_y= 0 + y \nabla^2 \zeta |_{y=0} + \frac{y^2}{2!} \frac{d}{dy} \nabla^2 \zeta |_{y=0} + ...$
$\displaystyle = 0 + 0 + 0+...$

But I'm struggling for a more formal proof... Could anyone point me in the right direction?