with(plots):with(LinearAlgebra):

EFG:=proc(X)
> local Xu,Xv,E,F,G;
> Xu:=<diff(X[1],u),diff(X[2],u),diff(X[3],u)>;
> Xv:=<diff(X[1],v),diff(X[2],v),diff(X[3],v)>;
> E:=DotProduct(Xu,Xu,conjugate=false);
> F:=DotProduct(Xu,Xv,conjugate=false);
> G:=DotProduct(Xv,Xv,conjugate=false);
> simplify([E,F,G]);
> end:

geoeq:=proc(X)
> local M,eq1,eq2;
> M:=EFG(X);
> eq1:=diff(u(t),t$2)+subs({u=u(t),v=v(t)},
> diff(M[1],u)/(2*M[1]))*diff(u(t),t)^2+subs({u=u(t),v=v(t)},
> diff(M[1],v)/(M[1]))*diff(u(t),t)*diff(v(t),t)-subs({u=u(t),
> v=v(t)},diff(M[3],u)/(2*M[1]))*diff(v(t),t)^2=0;
> eq2:=diff(v(t),t$2)-subs({u=u(t),v=v(t)},
> diff(M[1],v)/(2*M[3]))*diff(u(t),t)^2+subs({u=u(t),v=v(t)},
> diff(M[3],u)/(M[3]))*diff(u(t),t)*diff(v(t),t)-subs({u=u(t),
> v=v(t)},diff(M[3],v)/(2*M[3]))*diff(v(t),t)^2=0;
> eq1,eq2;
> end:

geoeq2:=proc(X)
> local M,G,E,F,D,Eu,Ev,Fu,Fv,Gu,Gv,eq1,eq2;
> M:=EFG(X);
> E:=M[1];F:=M[2];G:=M[3];
> D:=E*G-F*F;
> Eu:=diff(M[1],u);Ev:=diff(M[1],v);
> Fu:=diff(M[2],u);Fv:=diff(M[2],v);
> Gu:=diff(M[3],u);Gv:=diff(M[3],v);
> eq1:=diff(u(t),t$2)+subs({u=u(t),v=v(t)},
> (G*Eu-F*(2*Fu-Ev))/(2*D))*diff(u(t),t)^2
> +subs({u=u(t),v=v(t)},(G*Ev-F*Gu)/D)*diff(u(t),t)*
> diff(v(t),t)+subs({u=u(t),v=v(t)},
> (G*(2*Fv-Gu)-F*Gv)/(2*D))*diff(v(t),t)^2=0;
> eq2:=diff(v(t),t$2)+subs({u=u(t),v=v(t)},
> (E*(2*Fu-Ev)-F*Eu)/(2*D))*diff(u(t),t)^2+subs({u=u(t),
> v=v(t)},(E*Gu-F*Ev)/(D))*diff(u(t),t)*diff(v(t),t)
> +subs({u=u(t),v=v(t)},(E*Gv-F*(2*Fv-Gu))/(2*D))*
> diff(v(t),t)^2=0;
> eq1,eq2;
> end:

plotgeo:=proc(X,ustart,uend,vstart,vend,u0,v0,Du0,Dv0,
> T,N,gr,theta,phi)
> local sys,desys,u1,v1,listp,geo,plotX;
> sys:=geoeq(X);
> desys:=dsolve({sys,u(0)=u0,v(0)=v0,D(u)(0)=Du0,D(v)(0)=Dv0},
> {u(t),v(t)},type=numeric,output=listprocedure);
> u1:=subs(desys,u(t));v1:=subs(desys,v(t));
> geo:=tubeplot(convert(subs(u='u1'(t),v='v1'(t),X),list),
> t=0..T,radius=0.05,color=black,numpoints=N):
> plotX:=plot3d(X,u=ustart..uend,v=vstart..vend,grid=[gr[1],
> gr[2]],shading=XY,lightmodel=light3):
> display({geo,plotX},shading=XY,lightmodel=light2,
> scaling=constrained,orientation=[theta,phi]);
> end:

sphere:=<cos(u)*cos(v)|sin(u)*cos(v)|sin(v)>;

torus:=<(4+cos(u))*cos(v)|(4+cos(u))*sin(v)|sin(u)>;

EFG(sphere);

geoeq(sphere);

EFG(torus);

plotgeo(torus,0,2*Pi,0,2*Pi,Pi/2,0,0,1,30,150,[15,30],
> 105,33);