The error message is puzzling if you start integrating from t=0. It seems to indicate a misprint in the code, or a misinterpretation of how the package is supposed to be used.

1. I would first make sure that I can make the package work on a absolutely trivial problem (say two coupled linear first order equations, without any parameters).

2. Then I would start adding term by term to the equations, moving towards the equations of interest, and check that the solver continues to work in each case. It is generally a bad idea to just feed a complicated problem to a numerical package. How can you be sure that the solution is correct? How can you understand the properties of the solution?

3. Nonlinear initial value problems will usually develop singularities after some time. This is a property of the solution, which cannot be prevented by any solver package. So I would first try to integrate for a short time only.

Dear Kare

Thanks for your advise

yesterday with the help of Holy Spirit, finally I can find the errors and fix them.

enclosed is the final version.

Thanks again,

Yours,

Vic

Victor> ... with the help of Holy Spirit

I don't want to sound disrespectful or anything like that. But I have to tell the truth: Your solution is useless, and represents Mathematica in a pitiful way. How much do you pay per year for that crap? Why didn't the program warn you that the solution diverges after a time t=1.564269? After which the solution is just a demonstration of how good computers are at wasting time on 'garbage in--garbage out' processing.

It takes five minutes to code the solution in Python, which you can get for free (with a lot of freely available useful help on the net). A more time-consuming task is to make a decent figure of the solution. The best I can say of Mathematica is that they once hired a former master student of mine (but now he has found another job).

I tried to upload a figure of the solution here, but I only got error messages.

https://10.13140/RG.2.2.15146.82883

Code Find and plot solution of a Riccati equation

Mathematica is just a tool.

It is high time to try for an asymptotic  a = a_0/(t-t_)0) and b = b_0/(t-t_0). t_0 is arbitrary for the asymptotic and determined by initial conditions for the full problem. 

Solving the Riccati equations is just a beginning; these solutions must be inserted into new relations (shown in the papers by Ershkov) to obtain the velocity fields. Hence, diverging behaviour (simple pole singularities) of a(t), b(t)  does not necessarily imply unphysical behaviour. From a numerical point of view the problem becomes to locate each pole singularity on the t-axis, and get past it in a reliable way. That is a solvable technical problem. More important is to identify appropriate values for the free parameters in the general equation --- all set to 1 by  Victor (and me), and appropriate initial values for a and b. There is also a question of the dependence on positions (x, y, z), which has disappeared somewhere along the way.

I think this is an interesting alternative approach to the Navier-Stokes equations, although it is a bit unclear how far one may proceed analytically. And how much insight one can gain from the complicated analytic formulas. At some point the results have to be plotted to be understood; at that point it really does not matter how the numerical values are obtained. A calculation from some exceedingly complicated analytic expression is not necessarily any better than results found by solving the original Navier-Stokes equations numerically; it becomes a matter of numerical efficiency and reliability. But alternative numerical routes are always useful for checking results and methods.

Also, I don't think the Holy Spirit is the most trustworthy source for help on these matters...

You must specify the matrix type and probability distribution of matrix elements. The definition of sgn(Z) looks strange.

Try to check the prescription on a 1x1 (real) matrix.