#/***********************************************************
#orddif.rb -- 常微分方程式
#***********************************************************/

def sqr(x)  # $x^2$ 
return x * x
end

def F(x,  y)  # $F(x, y)$ 
return 1 - sqr(y)
end

def Fx(x,  y)  # $F_xend(x, y)$ 
return -2 * y * F(x, y)
end

def Fxx( x,  y)  # $F_xxend(x, y)$ 
return -2 * (sqr(F(x, y)) + y * Fx(x, y))
end

def euler( n,  nprint,  x0,  y0,  xn)  # Euler法 
x = x0;  y = y0;  h = (xn - x0) / n.to_f
for i in 1..n
    y += F(x, y) * h
    x = x0 + i * h
    if (i % nprint == 0); printf("% -14g % -14g\n", x, y); end
end
return y
end

def tayl3( n,  nprint,  x0,  y0,  xn)  # 3次Taylor級数 
x = x0;  y = y0;  h = (xn - x0) / n.to_f
for i in 1..n
    y += h * (F(x, y) + (h / 2) *\
                (Fx(x, y) + (h / 3) * Fxx(x, y)))
    x = x0 + i * h
    if (i % nprint == 0); printf("% -14g % -14g\n", x, y); end
end
return y
end

def runge4(n,  nprint,  x0,  y0,  xn)  # 4次Runge-Kutta法 
x = x0;  y = y0;  h = (xn - x0) / n.to_f;  h2 = h / 2.0
for i in 1..n
    f1 = h * F(x, y)
    f2 = h * F(x + h2, y + f1 / 2.0)
    f3 = h * F(x + h2, y + f2 / 2.0)
    f4 = h * F(x + h, y + f3)
    x = x0 + i * h
    y += (f1 + 2 * f2 + 2 * f3 + f4) / 6.0
    if (i % nprint == 0); printf("% -14g % -14g\n", x, y); end
end
return y
end

4.step(128, 2) do |n|
    printf("\nEuler法:          n = %d\n", n)
     euler(n, n / 4, 0, 0, 1)
    printf("\n3次Euler法:       n = %d\n", n)
     tayl3(n, n / 4, 0, 0, 1)
    printf("\n4次Runge-Kutta法: n = %d\n", n)
     runge4(n, n / 4, 0, 0, 1)
end
printf("\n正解\n")
for i in 1..4
    printf("% -14g % -14g\n", i / 4.0, Math::tanh(i / 4.0))
end

exit 0