> restart;

> alias(u=u(theta,r),v=v(theta,r),fr=fr(theta,r),ftheta=ftheta(theta,r),PHI=PHI(theta,r),h=h(hteta));

u, v, fr, ftheta, PHI, h

> alias(utheta=utheta(theta,r),ur=ur(theta,r),vtheta=vtheta(theta,r),vr=vr(theta,r));

u, v, fr, ftheta, PHI, h, utheta, ur, vtheta, vr

> EQN:=w*(vr+1/r^2*utheta+2/r*v+1/r^2*cos(theta)/sin(theta)*u+v*fr+1/r^2*u*ftheta)-(-1/r^3*v*u^2+v^2*vr+1/r^2*ur*u*v+1/r^2*u*v*vtheta+1/r^4*u^2*utheta)=0;

EQN := w*(vr+utheta/r^2+2*v/r+cos(theta)*u/(r^2*sin(theta))+v*fr+u*ftheta/r^2)+v*u^2/r^3-v^2*vr-ur*u*v/r^2-u*v*vtheta/r^2-u^2*utheta/r^4 = 0

> EQN:=subs(w=v^2+1/r^2*u^2,EQN);

EQN := (v^2+u^2/r^2)*(vr+utheta/r^2+2*v/r+cos(theta)*u/(r^2*sin(theta))+v*fr+u*ftheta/r^2)+v*u^2/r^3-v^2*vr-ur*u*v/r^2-u*v*vtheta/r^2-u^2*utheta/r^4 = 0EQN := (v^2+u^2/r^2)*(vr+utheta/r^2+2*v/r+cos(theta)*u/(r^2*sin(theta))+v*fr+u*ftheta/r^2)+v*u^2/r^3-v^2*vr-ur*u*v/r^2-u*v*vtheta/r^2-u^2*utheta/r^4 = 0

> u:=diff(PHI,theta):v:=diff(PHI,r):utheta:=diff(u,theta):ur:=diff(u,r):vtheta:=diff(v,theta):vr:=diff(v,r):

> PHI:=r-h(theta);

PHI := r-h(theta)

> EQN;

(diff(PHI, r)^2+diff(PHI, theta)^2/r^2)*(diff(PHI, `$`(r, 2))+diff(PHI, `$`(theta, 2))/r^2+2*diff(PHI, r)/r+cos(theta)*diff(PHI, theta)/(r^2*sin(theta))+diff(PHI, r)*fr+diff(PHI, theta)*ftheta/r^2)+di...(diff(PHI, r)^2+diff(PHI, theta)^2/r^2)*(diff(PHI, `$`(r, 2))+diff(PHI, `$`(theta, 2))/r^2+2*diff(PHI, r)/r+cos(theta)*diff(PHI, theta)/(r^2*sin(theta))+diff(PHI, r)*fr+diff(PHI, theta)*ftheta/r^2)+di...(diff(PHI, r)^2+diff(PHI, theta)^2/r^2)*(diff(PHI, `$`(r, 2))+diff(PHI, `$`(theta, 2))/r^2+2*diff(PHI, r)/r+cos(theta)*diff(PHI, theta)/(r^2*sin(theta))+diff(PHI, r)*fr+diff(PHI, theta)*ftheta/r^2)+di...

> simplify(%):

> %*r^4*sin(theta);

-r^2*diff(h(theta), `$`(theta, 2))*sin(theta)+2*r^3*sin(theta)-r^2*cos(theta)*diff(h(theta), theta)+r^4*fr*sin(theta)-r^2*diff(h(theta), theta)*ftheta*sin(theta)+3*diff(h(theta), theta)^2*r*sin(theta)...-r^2*diff(h(theta), `$`(theta, 2))*sin(theta)+2*r^3*sin(theta)-r^2*cos(theta)*diff(h(theta), theta)+r^4*fr*sin(theta)-r^2*diff(h(theta), theta)*ftheta*sin(theta)+3*diff(h(theta), theta)^2*r*sin(theta)...

> simplify(%);

-r^2*diff(h(theta), `$`(theta, 2))*sin(theta)+2*r^3*sin(theta)-r^2*cos(theta)*diff(h(theta), theta)+r^4*fr*sin(theta)-r^2*diff(h(theta), theta)*ftheta*sin(theta)+3*diff(h(theta), theta)^2*r*sin(theta)...-r^2*diff(h(theta), `$`(theta, 2))*sin(theta)+2*r^3*sin(theta)-r^2*cos(theta)*diff(h(theta), theta)+r^4*fr*sin(theta)-r^2*diff(h(theta), theta)*ftheta*sin(theta)+3*diff(h(theta), theta)^2*r*sin(theta)...

> subs(diff(h(theta),theta)=p(theta),%);

-r^2*diff(p(theta), theta)*sin(theta)+2*r^3*sin(theta)-r^2*cos(theta)*p(theta)+r^4*fr*sin(theta)-r^2*p(theta)*ftheta*sin(theta)+3*p(theta)^2*r*sin(theta)-p(theta)^3*cos(theta)+p(theta)^2*fr*r^2*sin(th...-r^2*diff(p(theta), theta)*sin(theta)+2*r^3*sin(theta)-r^2*cos(theta)*p(theta)+r^4*fr*sin(theta)-r^2*p(theta)*ftheta*sin(theta)+3*p(theta)^2*r*sin(theta)-p(theta)^3*cos(theta)+p(theta)^2*fr*r^2*sin(th...

> EQN2:=diff(p(theta),theta)=solve(%,diff(p(theta),theta));

EQN2 := diff(p(theta), theta) = -(-2*r^3*sin(theta)+r^2*cos(theta)*p(theta)-r^4*fr*sin(theta)+r^2*p(theta)*ftheta*sin(theta)-3*p(theta)^2*r*sin(theta)+p(theta)^3*cos(theta)-p(theta)^2*fr*r^2*sin(theta...EQN2 := diff(p(theta), theta) = -(-2*r^3*sin(theta)+r^2*cos(theta)*p(theta)-r^4*fr*sin(theta)+r^2*p(theta)*ftheta*sin(theta)-3*p(theta)^2*r*sin(theta)+p(theta)^3*cos(theta)-p(theta)^2*fr*r^2*sin(theta...

>