#This file finds the complex variable form of a real biharmonic function (Airy function) in the form f=g1+bar{g1}+z2*g2+z1*bar{g2} or f=2*Re{g1(z1)+z2*g2(z1)}.

#Replace the next line by the required real biharmonic (Airy) function in Cartesian coordinates, or comment it, specify the function phi(x,y) in the exterior part of the program and then call this file as a subroutine.


#phi:=x*y^3;

f:=expand(subs({x=(z1+z2)/2,y=I*(z2-z1)/2},phi));

#From (18.28), we have g_2^{\prime\prime} as

g2pp:=diff(f,z1,z1,z2);

#and hence, integrating,

g2p:=int(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:=int(g2p,z1);

#There is no need to introduce an additional constant here because any such term could be wrapped into g1.

g2b:=subs({z1=z2,I=-I},g2);

#In the next step, h=g1+g1b

h:=f-z2*g2-z1*g2b;

#Use the procedure rtoch for g1.

g1p:=diff(h,z1);

#A must be chosen to make g1p a function of z1 only

s:=coeff(g1p,z2);

sol1:=solve({s=0},{A});

g1p:=subs(sol1,g1p);
g2:=simplify(subs(sol1,g2));
g2b:=simplify(subs(sol1,g2b));

#Integrate to recover g1.

g1:=int(g1p,z1)+B;

g1b:=subs({z1=z2,I=-I},g1);

#Complex form of f from (18.27).

ff:=simplify(g1+g1b+z1*g2b+z2*g2);

err:=simplify(ff-f);

sol2:=solve({err=0},{B});

ff:=simplify(subs(sol2,ff));

test:=simplify(ff-f);

g1:=simplify(subs(sol2,g1));
g1b:=simplify(subs(sol2,g1b));
g2:=simplify(subs(sol2,g2));
g2b:=simplify(subs(sol2,g2b));