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_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?

Thanks for your help!!