The "pendulum" described above is sort of unphysical. The two huge masses can't stay at rest, they should be attracted by gravitation. So here is another system, three masses (restricted to a plane) but following Newtons law

n = 3;(* number of masses  *)
dim = 2;  (* dimension of space - here a plane  *)
pos = Table[x[i, j][t], {i, 1, n}, {j, 1, dim}];
mass = {10^-4, 1., 1.05};  (* Params of the system  *)
grav = 1.;  (* constant of gravity  *)

force[i_] := 
 - Sum[If[j == i,0, (grav mass[[j]])/ Sqrt[(pos[[i]] - pos[[j]]).(pos[[i]] - pos[[j]])] (pos[[i]] - pos[[j]])], {j, 1, n}]

Initial conditions

initpos = Thread /@ Thread[(pos /. t -> 0) == {{0, 5}, {-2, 0}, {2, 0}}];
initspeed = Thread /@ Thread[(D[pos, t] /. t -> 0) == {{0, .2}, {0, -.15}, {0, .15}}];  

ODE and its solution

sys = Flatten[Join[DGL, initpos, initspeed]];
sol = Flatten[NDSolve[sys, Flatten[pos /. a_[t] -> a], {t, 0, 200}]];
xx = Partition[(Flatten[pos]) /. sol, dim];

and look at it

Animate[
 Graphics[{PointSize[.05], Point /@ xx /. t -> tt}, 
  PlotRange -> {{-6, 6}, {-10, 10}}],
 {tt, 0, 50}]

It seems that the small mass finally picks up so much speed that it leaves the system

I thought about your problem and - may be I am wrong - would say

If you have two huge masses M positioned at { - x, 0 } and { x, 0 } and if these masses are effectifely immobile and if you have a tiny mass m at { 0, h } these three form an isoceles triangle and I get as equation of motion for m

h''[ t ]  == -  G M h[ t ] /  Sqrt[ x^2 + h[ t ]^2 ]^3

So, where is your 2nd equation or what do you mean by it ?

The system described before is essentially sort of a pendulum: the tiny mass oscillates between the two huge masses:

Clear[h]
sol = NDSolve[{h''[t] == -h[t]/Sqrt[1 + h[t]^2]^3, h[0] == 1,  h'[0] == 0}, h, {t, 0, 10}]
h = h /. Flatten[sol]
Plot[h[t], {t, 0, 10}]

the first problem (for me) is the y''. i know it's possible to reduce an any-order system to a first level system but do not explain or show proof here.

substitution may be helpful

if i let

ksi=1

u==(1/(y^2 + 1/4))^(3/2)

y' == -1+1/2u

y'' == (1 - 1/u) y

the first diffeq is a linear non-homogenous first order ordinary differential equation linear in u, F(u,y')=0 and F(u,y'')=Q(x) : when solved u will have to be substituted and resolved.

the second is a second order equation of the same properties but it is homogenous (and needs u substituted afterward).

i haven't tried solving either of the above: they may be too complicated to solve by hand or non-linear after u is found, idk

once either is solved it could be substituted for in the other equation by a grab bag of means. as y (with precaution), finding y' one can find from it y'' and substitute that. doing so might make the system easier to understand

that's as far as i go - i got homework to do :)

i can say if you solve a system of linear equations using series method then you'd prefer having initial conditions that method needs, where 0 suggests maclaurin series, and that ultimately NDsolve and numeric approximation should be understood because it can do allot for you that DSolve cannot: if approximation methods understood, if not NDSolve might serve to confuse.

i have not reviewed the ODE properties of the 2 equations i'm not swift at reading them in mathematica yet

i suggest if there is some discussion over validity that the original poster plots (isoclines) of the two (differential equations) and see visually (though it's possible to be "tricked" by such plots)

https://reference.wolfram.com/language/howto/PlotTheResultsOfNDSolve.html

EquationTrekker

i would say if Hans Dolhaine said there is a problem then to believe him

As I consider the restricted three-body problem, three bodies form an equilateral triangle. Variable h(t) is height of this equilateral triangle and massless body is on vertex of this triangle. The differential equations system describes the dynamics of the restricted three-body problem.

The system is inconsistent and correct. I expect to plot numerical values of functions h[t] and u[t].

Thanks a lot again. Actually I am doing PhD in celestial mechanics, and this is dynamical equations of the restricted three-body problem. I am From Almaty, Kazakhstan.

I need to solve the differential equations system and obtain plots of h[t] and u[t].

Trying to work out :)

But you can do the following. Your system is incosistent, but you could ask for another function u. Redefine your system and solve it numerically (NDSolve)

ksi = 1;
system = {
   h'[t] == ksi/2*((1/(h[t]*h[t] + 1/4)^3/2) - 1),
   h[0] == 0,
   u'[t] == -((1/(h[t]*h[t] + 1/4)^3/2) - 1)*h[t],
   u[0] == 0
   };
sol = NDSolve[system, {h, u}, {t, 0, 10}] // Flatten;

Get the functions h and u

h = h /. sol;
u = u /. sol;

And their derivatives

hp = Derivative[1][h ];
up = Derivative[1][u];

Calculate values of these at certain values of t

hppts = Table[{t, hp[t]}, {t, 0, 2, .05}];
uppts = Table[{t, up[t]}, {t, 0, .6, .05}];

and look whether these values fit to the right-hand sides of your differential equation(s)

Plot[ksi/2*((1/(h[t]*h[t] + 1/4)^3/2) - 1), {t, 0, 1}, 
 PlotRange ->All,  Epilog -> {Red, PointSize[.015], Point /@ hppts}]

and

Plot[-((1/(h[t]*h[t] + 1/4)^3/2) - 1)*h[t], {t, 0, .6}, 
 PlotRange -> All, Epilog -> {Red, PointSize[.015], Point /@ uppts}]

Seems your system is incosistent:

With

hprime = 1/2 (-1 + 1/(2 (1/4 + h[t]^2)^3));
hdoubleprime = h[t] (1 - 1/(2 (1/4 + h[t]^2)^3));

and

FullSimplify[Together[D[hprime, t] - hdoubleprime]]

I don't get zero

I will suggest a few things.

(1) Clear all Global context symbols (or restart the kernel).

(2) Remove the second derivative equation. It makes the system overdetermined.

(3) If you want to plot a solution, give an initial value in the system.

Beyond that, it is possible the system has no (known) exact solution, so be prepared to use NDSolve instead.

What have you tried? What is unclear about the error message?