(*                     urθ                     *)

(*This program calculates the stresses and displacements for a given biharmonic stress function ϕ in Cartesian coordinates. To test it, uncomment the next line.*)

ϕ=A*r^2*Log[r]

(*This will calculate the displacements for plane stress.For plane strain,uncomment κ=3-4*nu and comment κ=(3-nu)/(1+nu). Alternatively,you can comment both lines and get the resuts in terms of κ.*)
(*κ=3-4*ν*)
κ=(3-ν)/(1+ν)

σrr=Simplify[(1/r)*D[ϕ,r]+(1/r^2)*D[ϕ,θ,θ]]
σθθ=D[ϕ,r,r]
σrθ=-Simplify[D[D[ϕ,θ]/r,r]]

err=Expand[((κ+1)*σrr-(3-κ)*σθθ)/4]
eθθ=Expand[((κ+1)*σθθ-(3-κ)*σrr)/4]
erθ=Expand[σrθ]

tmur=Integrate[err,r]
e1=Expand[r*eθθ-tmur]
tmuθ=Expand[Simplify[Integrate[e1,θ]]]

tmur1=tmur+D[f[θ],θ]
tmuθ1=tmuθ+g[r]-f[θ]

h=Expand[Simplify[r*(erθ-(1/2)*((1/r)*D[tmur1,θ]+D[tmuθ1,r]-tmuθ1/r))]]

(*The following two lines are a trick to separate the functions of r from the functions of θ.*)

h1=Integrate[D[h,θ],θ]
h2=Simplify[h-h1]

DSolve[h1==0,f[θ],θ,GeneratedParameters->A]
f1[θ]=f[θ]/.%
tmur2=tmur+D[f1[θ],θ]
tmuθ2=tmuθ+g[r]-f1[θ]
DSolve[h2==0,g[r],r,GeneratedParameters->B]
uθ1=tmuθ2/(2*μ)/.%
ur1=tmur2/(2*μ)

ur=Expand[ur1]
uθ=Expand[Simplify[uθ1]]