Cauchy Reimann equations-help with formal proof..

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

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

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

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

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

$p_x = p_x|_{y=0} +y p_{xy}|_{y=0} + \frac{y^2}{2!} p_{xyy}|_{y=0} + ...$
$\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 $y=0$,

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

and so on..

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

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