What happens if you replace all your decimal approximations with exact fractions so there will be no round-off error and you write your own simple Euler numerical solver? Your equations give explicit expressions for your second derivatives. It looks like you can get linear approximations of your first derivatives from those and then get linear approximations of your functions from those. Iterate that over each small step and see how it behaves.

I understand many users want a magic "do what I want" button. I realize there must be a better way, but if you can't understand and successfully control the inner works of any of the black box methods then writing a very simple numerical solver that you could understand might get you an answer today and that might be enough to get you started. That only seems feasible because the form of your system and initial conditions are simple enough to perhaps allow this. I've written a dozen lines to try to do that but I'm still trying to convince myself it is correct.

Dear Bill. Yes I was able to make work after I restricted the time to very small number. I tried your suggestion and it work too. Thank you. The only problem causing the code not work is the nonlinear power of 3/2 in the equations F12 and F23. If you set 3/2 to 1 "cancel the nonlinearity" in the system everything works perfectly but again it's not going to be able to give the results I am looking for since I linearized my system. I am not sure if there is any other way of solving since I have been trying to for a few days. Thank you again.