Almost but no.

The other equations are

2sin(theta) phi' R' + 2R cos(theta) phi' theta ' + Rsin(theta) phi '' = 0

2R' theta ' +R theta'' -Rsin(theta) cos(theta) phi' ^2 = 0

is this too hard to do analytically?

The equations are actually

x'' + x e^(-R^2)=0

z'' + z e^(-R^2)=0

y'' + y e^(-R^2)=0

where x'' and y'' and z'' are in equations 75 in :

http://mathworld.wolfram.com/SphericalCoordinates.html

How do I do this now?

for R theta and phi

> with(plots);
> graph1 := plot([sqrt(4*r^6-2*r^4), 3*r^4-2*r^2, r = 1/2 .. 10], -2 .. 2, -.5 .. 5, labels = [j, h]); graph2 := plot([-sqrt(4*r^6-2*r^4), 3*r^4-2*r^2, r = 1/2 .. 10], -2 .. 2, -.5 .. 5, labels = [j, h]); origin := [[0, 0]];
                                  [[0, 0]]
> ptplot := plot(origin, style = POINT, color = BLACK, symbol = BOX);
> display(graph1, graph2, ptplot);
 

Nevermind.

Is this ever true?

Sorry, how do I assume 'a' is nonzero and solve for the inequality to avoid the error

How do I assume a is nonzero?

won't work, d you know why?

`assuming`([solve({b^2/a^2-2*b/a+2*b+1-2*a+a^2 > 0, a-b < 0}, [a, b])], [a, 'real', a < 0, b, 'real', 0 < b])

haha

And how would I draw multiple solution curves along with the vector field?

What I had written in the OP didn't resemble what i was expecting i.e. a hyperbola in addition to the null cline that was given in your command, hyperb