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, 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^2;
   Eu:= diff(E, u);   Ev:= diff(E, v);
   Fu:= diff(F, u);   Fv:= diff(F, v);
   Gu:= diff(G, u);   Gv:= diff(G, 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: