As is often necessary on a problem like this, it helps to look carefully at what you have.
x = sym('x',[1 2]);
mu_earth = 398600.441;
R_earth = 6378;
M_dot_pert = 1.11337451*10^(-3);
Omega_dot_pert = 1.991063853*10^(-7);
J2=0.0010826269;
F(1) = Omega_dot_pert +1.5*J2*((R_earth/x(1))^2)*M_dot_pert*cos(x(2))
F(2) = M_dot_pert - sqrt((mu_earth/x(1)^3))*(1+1.5*J2*(((R_earth/x(1))^2)*(1-1.5*(sin(x(2))^2))))
So you have some REALLY HUGE numbers in there. And you are raising large numbers (on the order of 7000, to powers. That alone MAY be a significant problem.
As bad is the idea that you are forming the sin and cos of numbers that are on the order of 100.
Is the number 97 in DEGREES????????? If it is in degrees, then you clearly do not understand that sin and cos assume RADIANs. If you want to work in degrees, then you need to use sind and cosd. In fact, I can be almost positive you think that 97 is a number in degrees, since that is perhaps the most common mistake we ever see on this site.
97 degrees is
deg2rad(97.331)
So 1.6987 RADIANS. Again, if you want to work in degrees, you will need to write the above expressions using sind and cosd. Take your pick. I've done there here. Note that the symbolic toolbox converts degrees to radians, then uses sin and cos.
F(1) = Omega_dot_pert +1.5*J2*((R_earth/x(1))^2)*M_dot_pert*cosd(x(2))
F(2) = M_dot_pert - sqrt((mu_earth/x(1)^3))*(1+1.5*J2*(((R_earth/x(1))^2)*(1-1.5*(sind(x(2))^2))))
But at the initial point, now we can look at your objective.
x0 = [6858,97.331];
vpa(subs(F,x,x0),5)
And we see here that it results in already very small numbers, near the default tolerance for fsolve.
If I compute the gradient, I'd bet that again, we will see small numbers. You will need to learn to scale your problems, using units that are appropriate, yet will not cause numerical problems. Effectively, you need to learn why it is that astronomers work in light years and mega-parsecs, not in feet or meters.
