(* This file finds the complex variable form of a real harmonic function, f in the form f=g+bar{g} or f=2*Re{g} (see equation (18.24)). *)
(* Replace the next line by the required real harmonic  function in Cartesian coordinates, or delete it,specify the function f(x,y) in the exterior part of the program and then call this file as a subroutine. *)
f=x*y
f=f/.{x->(z1+z2)/2,y->I*(z2-z1)/2}

(* Evaluate g^prime from equation (18.26).*)
gp=D[f,z1]
(* Integrate to recover g *)
g=Integrate[gp,z1]+A
(* Construct the conjugate of g *)
gb=g/.{z1->z2,z2->z1,I->-I}
(* Complex form of f from (18.24). *)
ff=Simplify[g+gb]
err=Simplify[ff-f]
sol=Solve[{err==0},{A}]
g=g/.sol
gb=gb/.sol