гравітаційне прискорення

Ви прораховували поведінку апарату при русі навколо Місяця. Проте поведінка Місяця ви не прораховували. От і все.
Взагалі-то, маса Місяця настільки велика, що і ніхто її НЕ прораховує, а очманілі метеорити обурюють поведінку Місяця більше, ніж нещасний кораблик.

ще раз і повільно. При русі по витягнутому еліпсу швидкість в нижній точці еліпса вище, ніж в бічних або у верхній. Це нормально, апарат при цьому ні на які параболи / гіперболи не переходить. Наближаємося до планети - швидкість зростає, віддаляємося - падає назад. Все строго симетрично щодо довгої осі еліпса.

І ще питання - у вас маса апарату порівнянна з масою планети або незначна, як в реальному житті. Щодо першого випадку є одна підозра.

Теж саме - ще раз і повільно.
Маси Землі і Місяця брав з довідника.
Маса корабля - 10 тонн.
Корабель встановлювався в точку на висоті 300 000 км від центру землі і йому повідомлялася тангенціальна швидкість розміром 1020 м / с.
Корабель виходив на практично кругову орбіту і робив там десятка півтора оборотів.
При кожному оберті він проходив повз Місяця яка трохи змінювала його орбіту своїм тяжінням.
Зрештою орбіта його руйнувалася і він полетів у бік Землі і врізався в неї.
Зробив іншу швидкість - 1019 м, с.
Рух корабля було практично таке ж, але на цей раз він не врізався в Землю, а пролітав повз неї, отримуючи потужний гравітаційне прискорення

Якщо у вас місяць і земля нерухомі - значить насправді у вас все відбувалося не в інерціальній системі відліку а в системі обертається з частотою 1об / 4неделі. Вже одне це змушує засумніватися у вірності математики.

Місяць прибрати не пробували? Що в цьому випадку відбувається?

Здається, я здогадався про джерело прискорення в програмі, що використовує чисельні методи.
Припустимо, тіло виявилося в точці т0, яка знаходиться в максимальній близькості до центру планети.
У цій точці у тіла мінімальна потенційна і максимальна кінетична енергія, отже, максимальна швидкість.
Поки все в порядку - все по законам небесної механіки.
Далі програма починає обчислювати наступну точку т1.
Вона, очевидно, знаходиться на більшій відстані від центру планети, так як т0 було на мінімальній відстані.
І в цій точці, за законами небесної механіки, швидкість повинна бути МЕНШЕ.
Але програма всуває в точку т1 швидкість точки т0, що, строго кажучи, незаконно.
І швидкість в цій точці виходить БІЛЬШЕ, ніж повинна бути за законами небесної механіки.
Ось і джерело збільшення швидкості.
У разі, коли тіло далеко від центру планети, ця помилка занадто мала, щоб бути помітною.
А ось поблизу центру тяжіння коли і сили і швидкості від точки до точки, різко, іноді на порядок, змінюються, вона і вилазить.
Так що вважаю, що тема закрита.

Думаю, ви маєте рацію. Істотна помилка розрахунку проявляється имено поблизу масивних тіл, коли на кожному кроці розрахунку істотно змінюється потенціал гравітації.
Я використовував метод, що дає дещо меншу похибку. Розраховував зміна гравітаційного потенціалу при переміщенні тіла в передбачувану точку Z1 за час DT, в кінці періоду + DT пересчітавал зміна імпульсу тіла, що направляє вектора імпульсу обчислювалася як рівнодіюча сил тяжіння мас системи в точці між Z0 і Z1. Приблизно так.

використовуйте чисельну схему v1 = v0 + 0.5 * dt * (a0 + a1); r1 = r0 + 0.5dt * (v0 + v1), де v - швидкості, a - прискорення, r - радіус-вектор, 0 - величини в початковий момент часу, а 1-- в новий. на першій ітерації прирівняти 1е до 0м, а потім пару раз виженете, за декілька ітерацій рішення зійдеться. ну і так на кожному моменті часу.