(* This file converts a real Airy stress function phi(x,y) into the two complex potentials chi(z1),theta(z1).In this file we write chi(z1)=chi and its conjugate as chib,etc. *)
ϕ=x^2
f=Expand[ϕ/.{x->(z1+z2)/2,y->I*(z2-z1)/2}]

(* From (18.28),we have g_2^{\prime\prime} as *)
g2pp=D[f,z1,z1,z2]
(* and hence, integrating twice *)
g2p=Integrate[g2pp,z1]+A
(* The constant in g2p+g2pb must be real,so there is no loss in generality in taking A to be real.Integrate again to get *)
g2=Integrate[g2p,z1]
(* There is no need to introduce an additional constant here because any such term could be wrapped into g1. *)
(* The conjugate is *)
g2b=Simplify[Conjugate[g2],Assumptions->{Element[z1,Reals],Element[A,Reals]}]/.z1->z2

(* In the next step,h=g1+g1b *)
h=f-z2*g2-z1*g2b
(* Use the procedure rtoch for g1. *)
g1p=D[h,z1]

(* A must be chosen to make g1p a function of z1 only. *)
s=Coefficient[g1p,z2]
sol=Solve[{s==0},{A}]
g1p=g1p/.sol
g2=g2/.sol
g2b=g2b/.sol

(* Integrate to recover g1. *)
g1=Integrate[g1p,z1]
g1b=Simplify[Conjugate[g1],Assumptions->{Element[z1,Reals]}]/.z1->z2

ff=Simplify[g1+g1b+z1*g2b+z2*g2]

(* Here we write chi(z1)=chi,theta(z1)=theta and the conjugates as chib,thetab *)
χ=2*g2
χb=2*g2b
θ=2*D[g1,z1]
θb=2*D[g1b,z2]
(* We can also then calculate the displacements as *)
tmu=κ*χ-z1*D[χb,z2]-θb