# dy/dx=1/(x+sqrt(x^2+1))-y/sqrt(x^2+1)
# 2nd-order RK method
import numpy as np
import matplotlib.pyplot as plt
def f(y, x):
return 1/(x+np.sqrt(x**2+1))-y/np.sqrt(x**2+1)
def RK2(f, x, t):
h=t[1]-t[0]
for i in range(len(t)-1):
k_1=h*f(x[i], t[i])
k_2=h*f(x[i]+k_1/2, t[i]+h/2)
x[i+1]=x[i]+k_2
return x
if __name__=="__main__":
x_i=0; y_i=0
x_f=1
x=np.linspace(x_i, x_f, 1000)
y=np.zeros(len(x), float)
y[0]=y_i
y=RK2(f, y, x)
plt.plot(x, y)
plt.show()