Метод Рунге-Кутты для численного решения нелинейного дифференциального уравнения второго порядка - Octave

Итак, у меня есть это аналитически неразрешимое нелинейное дифференциальное уравнение второго порядка:

x '' = (- 1 / x) ((x ') ² + gx-gh + ax')

к которому я некоторое время пытался применить метод Рунге-Кутта в Octave, но не получил хороших результатов. Это код, который я использовал:

graphics_toolkit('gnuplot')

to=0
tf=2.8
N=150
dt=[tf-to]/N

b=2.65
g=9.81
h=0.07

y(1)=0.015
z(1)=0
t(1)=to
for i=1:N-1
  y(i+1)= y(i) - dt.*z(i)
  z(i+1)= z(i) - dt.*[z(i).**2 + g.*y(i) - g.*h + b.*z(i)]./y(i)
  t(i+1)= t(i) + dt
endfor

plot(y,t)
title=('Numerical solution')
xlabel=('Time')
ylabel=('Position')

Честно говоря, я не могу сказать, в чем проблема, может быть, я что-то не понимаю в методе, или, может быть, в моем коде отсутствует что-то важное.

Спасибо за ваше время и внимание.