# dy/dx=1/(x+sqrt(x^2+1))-y/sqrt(x^2+1)
# 4th-order RKM HW14
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 RK4(f, x, t):
h=t[1]-t[0]
for i in range(len(x)-1):
c_1=h*f(x[i], t[i])
c_2=h*f(x[i]+c_1/2, t[i]+h/2)
c_3=h*f(x[i]+c_2/2, y[i]+h/2)
c_4=h*f(x[i]+c_3, t[i]+h)
x[i+1]=x[i]+(c_1+2*c_2+2*c_3+c_4)/6
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_like(x, float)
y[0]=y_i
y=RK4(f, y, x)
plt.plot(x, y)
plt.show()