syms phi theta r f(r,theta,phi)

Df(1)=diff(f,r);
Df(2)=diff(f,theta);
Df(3)=diff(f,phi);


E(1,1)=cos(theta)*sin(phi); E(2,1)=sin(theta)*sin(phi); E(3,1)=cos(phi); 

E(1,2)=-sin(theta); E(2,2)=cos(theta); E(3,2)=0; 

E(1,3)=cos(theta)*cos(phi); E(2,3)=sin(theta)*cos(phi); E(3,3)=-sin(phi); 


J(1,1)=cos(theta)*sin(phi); J(1,2)=sin(theta)*sin(phi); J(1,3)=cos(phi); 

J(2,1)=-r*sin(phi)*sin(theta); J(2,2)=r*cos(theta)*sin(phi); J(2,3)=0; 

J(3,1)=r*cos(theta)*cos(phi); J(3,2)=r*sin(theta)*cos(phi); J(3,3)=-r*sin(phi); 

simplify(inv(E)*inv(J)*Df.')