Hello, I'm having trouble getting NDSolve to solve a system of DEs. Here's my code below, can you give me an idea of what I need to correct please? THANKS SO MUCH!!!

*User inputs chart functions. ee = x (u, v); ff = y (u, v);

 gg = cos (v);*)

ee = Cos[u]*Sin[v]
ff = Sin[u]*Sin[v]
gg = Cos[v]

eeh[u_, v_] := ee[x][u, v];

ffh[u_, v_] := ff[x][u, v];

ggh[u_, v_] := gg[x][u, v]


(*User inputs min and max values for u and v*)

u0 = -Pi
u1 = Pi

v0 = -Pi
v1 = Pi

(*(User inputs  initial p0 and q0)*)

p0 = 1
q0 = 1

tmin = 0
tmax = 200

a = mat = {{D[ee, u], D[ee, v]}, {D[ff, u], D[ff, v]}, {D[gg, u], 
    D[gg, v]}}

(a) // MatrixForm
Print["Derivative Matrix"]

e = D[ee, u]*D[ee, u] + D[ff, u]*D[ff, u] + D[gg, u]*D[gg, u]
e = TrigReduce[e]

h = D[ee, u]*D[ee, v] + D[ff, u]*D[ff, v] + D[gg, u]*D[gg, v]
h = TrigReduce[h]

g = D[ee, v]*D[ee, v] + D[ff, v]*D[ff, v] + D[gg, v]*D[gg, v]
g = TrigReduce[g]

bb = mat = {{e, h}, {h, g}}
bb // MatrixForm
Print["Metric, g"]
bbb = Inverse[bb]
Print["Christoffel Equations: gamma_i,j,k"]

b = 0.5*(Part[bbb, 1, 1]*D[e, u] + 
    Part[bbb, 1, 2]*(2*D[h, u] - D[e, v]))
Print["gamma_1,1,1"]

c = 0.5*(Part[bbb, 2, 1]*D[e, u] + 
    Part[bbb, 2, 2]*(2*D[h, u] - D[e, v]))

Print["gamma_1,1,2"]

d = 0.5*(Part[bbb, 1, 1]*D[e, v] + Part[bbb, 1, 2]*D[g, u])

Print["gamma_1,2,1"]

i = 0.5*(Part[bbb, 2, 1]*D[e, v] + Part[bbb, 2, 2]*D[g, u])
Print["gamma_1,2,2"]

j = 0.5*(Part[bbb, 1, 1]*D[e, v] + Part[bbb, 1, 2]*D[g, u])
Print["gamma_2,1,1"]

k = 0.5*(Part[bbb, 1, 1]*(2*D[h, v] - D[g, u]) + 
    Part[bbb, 1, 2]*D[g, v])
Print["gamma_2,2,1"]

l = 0.5*(Part[bbb, 2, 1]*D[e, v] + Part[bbb, 2, 2]*D[g, u])
Print["gamma_2,1,2"]

m = 0.5*(Part[bbb, 2, 1]*(2*D[h, v] - D[g, u]) + 
    Part[bbb, 2, 2]*D[g, v])
Print["gamma_2,2,2"]

(* gdesolv[u0,v0, p0, q0, t0, tmin, tmax]:=*)
soln4 = NDSolve[ {u'[t] == p[t], v'[t] == q[t],
   p'[t] + b*p[t]^2 + 2*c*p[t]*q[t] + i*q[t]^2 == 0,
   q'[t] + j*p[t]^2 + 2*l*p[t]*q[t] + m*q[t]^2 == 0, u[t0] == u0, 
   v[t0] == v0, p[t0] == du0, q[t0] == dv0},
  {u[t], v[t], p[t], q[t]}, {t, tmin, tmax}]

ParametricPlot[Evaluate[{u[t], v[t]} /. soln4], {t, tmin, tmax}]



ParametricPlot3D[{ee, ff, gg}, {u, u0, u1}, {v, v0, v1}]