多変数の4次のルンゲクッタのプログラムを書いているのですが、どうもうまくいきません。

今、厳密解を

y(1) = exp(x) + exp(-x)
y(2) = exp(x) - exp(-x)

とおいて、その微分方程式

dydx(1) = exp(x) - exp(-x)
dydx(2) = exp(x) + exp(-x)

を数値的に解いて、厳密解と比較しています。

サブルーチンのderive中で、

dydx(1) = exp(x) - exp(-x)
dydx(2) = exp(x) + exp(-x)

とおいた場合は、きちんと数値解と厳密解が一致するのですが、

dydx(1) = y(2)
dydx(2) = y(1)

とおくと、両者は一致しません。これはなぜでしょうか？
考えているのですが、一向にわかりません。
プログラムは以下になります。よろしくお願いします。