(*  CONICAL BAR LOADED BY A MOMENT ABOUT THE x-AXIS  *)(*  Stresses vary with R^{-3} so stress functions must vary with R^{-1} (A and E) and R^{-2} (B).The required function for B is R^{-2}P_1^1(cos(beta))sin(theta),but this form would give m>n for A,E,so we need to use the log(R+z) form of equations (22.37).To get the appropriate non-axisymmetric function,take log(R+z) and differentiate with respect to y,obtaining  *)

ϕ=C[1]*Sin[β]*Sin[θ]/(R*(1+Cos[β]))
ω=C[2]*Sin[β]*Sin[θ]/R^2
ψ=C[3]*Sin[β]*Cos[θ]/(R*(1+Cos[β]))

σRR=Simplify[D[ϕ,R,R]+R*Cos[β]*D[ω,R,R]-2*(1-ν)*Cos[β]*D[ω,R]+(2*ν/R)*Sin[β]*D[ω,β]+(2/R)*D[ψ,R,θ]-(2/R^2)*D[ψ,θ]]
σθθ=Simplify[(1/R)*D[ϕ,R]+(Cot[β]/R^2)*D[ϕ,β]+(1/(R*Sin[β])^2)*D[ϕ,θ,θ]+(1-2*ν)*Cos[β]*D[ω,R]+(2*ν/R)*Sin[β]*D[ω,β]+(Cos[β]^2/(R*Sin[β]))*D[ω,β]+(Cot[β]/(R*Sin[β]))*D[ω,θ,θ]+(2/(R*Sin[β])^2)*D[ψ,θ]-(2/R)*D[ψ,R,θ]-(2*Cot[β]/R^2)*D[ψ,β,θ]]
σββ=Simplify[(1/R)*D[ϕ,R]+(1/R^2)*D[ϕ,β,β]+(Cos[β]/R)*D[ω,β,β]+(1-2*ν)*Cos[β]*D[ω,R]+(2*(1-ν)*Sin[β]/R)*D[ω,β]+(2*Cot[β]/R^2)*D[ψ,θ,β]-2*(Cot[β]/R)^2*D[ψ,θ]]
σθβ=(1/(R^2*Sin[β]))*D[ϕ,θ,β]-(Cot[β]/(R^2*Sin[β]))*D[ϕ,θ]+(Cot[β]/R)*D[ω,β,θ]+(2*(1-ν)/R)*D[ω,θ]-(1/(R*Sin[β]^2))*D[ω,θ]+(Cot[β]/(R^2*Sin[β]))*D[ψ,θ,θ]-(Sin[β]/R)*D[ψ,R,β]-(Cos[β]/R^2)*D[ψ,β,β]+(1/(R^2*Sin[β]))*D[ψ,β]
σβR=(1/R)*D[ϕ,β,R]-(1/R^2)*D[ϕ,β]+(1-2*ν)*Sin[β]*D[ω,R]+Cos[β]*D[ω,β,R]-(2*(1-ν)*Cos[β]/R)*D[ω,β]+(1/R^2)*D[ψ,θ,β]+(Cot[β]/R)*D[ψ,θ,R]-(2*Cot[β]/R^2)*D[ψ,θ]
σRθ=(1/(R*Sin[β]))*D[ϕ,R,θ]-(1/(R^2*Sin[β]))*D[ϕ,θ]+Cot[β]*D[ω,R,θ]-(2*(1-ν)*Cot[β]/R)*D[ω,θ]+(1/(R^2*Sin[β]))*D[ψ,θ,θ]-Sin[β]*D[ψ,R,R]-(Cos[β]/R)*D[ψ,R,β]+(2*Cos[β]/R^2)*D[ψ,β]+(Sin[β]/R)*D[ψ,R]

t1=σβR/.β->β0
t2=σββ/.β->β0
t3=Expand[σθβ/.β->β0]

c1=Simplify[Expand[t1/Sin[θ]]]
c2=Simplify[t2/Sin[θ]]
c3=Simplify[t3/Cos[θ]]

t=Simplify[σβR*Sin[θ]+σRθ*Cos[θ]*Cos[β]]
c4=Integrate[t, {θ, 0,2*π}]
c5=Simplify[Integrate[c4*R^3*Sin[β], {β, 0,β0}]]

Simplify[Solve[{c1==0,c2==0,c3==0,c5==M},{C[1],C[2],C[3]}]]

σRR=Simplify[σRR/.%20]

σθθ=Simplify[σθθ/.%20]

σββ=Simplify[σββ/.%20]

σθβ=Simplify[σθβ/.%20]

σβR=Simplify[σβR/.%20]

σRθ=Simplify[σRθ/.%20]