please help me by finding the bug to make the program stable. <A HREF="http://www.problem.cfd">www.problem.cfd</a>

the correct page is: http://www.wh12.tu-dresden.de/~andre

Probably it is code bug (I have not looked at your program - no comments and variable names in German ) but more similar that you have met problem with calculation stability problem in critical section of Laval nozzle. As it is described in many textbooks (Pirumov-Roslyakov, Sternin,...) the steady subsonic-supersonic quasi-one-dimensional Euler system of equations can be rewritten as one equation for M and at critical cross-section in case of frozen flow (or after critical cross-section in case of non-equilibrium flow) it is special point with term DF/dx /(M**2-1). Some solvers (like LSODE) allow to pass this point smoothly, othes - not. One of possible methods to avoid this problem - to use "inverse" method - when you set not F(x), but ro(x), which can be calculated from isoentrophic relations (even in case of non-equilibrium flows density distribution is very close to frozen distribution). In this case this special point dissapear.

(1). Most density-based methods require some kind of artificial viscosity terms to reduce the oscillations in the solution. (2). My suggestion is: pick a well-known technical paper from the Journal, and try to follow the working method. (3). Book by Dale Anderson, John Tannehill and Richard Pletcher, "computational fluid mechanics and heat transfer" is a good starting point.

You may find it useful to reference the CFD book by John D Anderson. He discusses this problem in detail