Hi Luis,
I am going to assume that r is constant and not a function of theta. One way that could happen is if these equations are modeling small excursions about a fixed r, is that the case?
Anyway, if r is constant, then then an analytic solution is possible. Here is an outline, without doing out the details. Lots of small steps to go through, which is typical. First, in order for the units to work out it looks like both b and mu have units of momentum^2. After multiplying the first equation by (1/h)*r^(-1/2) and the second by (mu/h)*r^(-3/2) and rearranging, you can arrive at
c1*ga' +c2*ga''' = c3*gb +c4*gb''
c5*gb' +c6*gb''' = c7*ga +c8*ga''
where the c's are uncomplicated factors that you can work out. Some of the c's include two dimensionless constants that set the problem,
C = l/(h*r^2) and D = mu/h.
Since all terms include a linear factor of ga or gb, their size is undetermined within an overall multiplicative constant, which could have dimensions. Without including that constant, an analytic solution is
ga = cos(m*theta+phi) gb = A*sin(m*theta+phi) (to solve for m and A)
Phi can be determined later by initial conditions (maybe the overall size can be too). After plugging this in and doing the derivatives of the trig functions, you will get a couple of equations
-c1*m +c2*m^3 = A*(c3 -c4*m^2)
A*(c5*m -c6*m^3) = c7 -c8*m^2
Multiplying the two equations together to eliminate A gives a cubic polynomial in the variable m^2. At this point if not before, you can bring in Matlab to solve this. There will be six possible values for m. For each one you can back substitute into one of the equations to find the appropriate A.
Depending on the constants it could be that m^2 has one positive root which would allow two related physical solutions m = +-sqrt(root) for ga and gb.
Best Answer