I have to solve numerically the following initial value problem, which physical origin is discussed in this article by G.I. Taylor (1941). Here $\eta$ is the independent variable and ranges beetween $0$ and $1$, $\gamma$ is fixed and $=1.4$ in my simulation: $$f'=\dfrac{f[-3\eta + (3+\gamma/2)\phi - 2\gamma \dfrac{\phi ^2}{\eta}]}{(\eta-\phi)^2-f/\psi},$$$$\phi'=\dfrac{\frac{1}{\gamma} f'}{\psi(\eta-\phi)}-\dfrac{3}{2}\dfrac{\phi}{\eta-\phi},$$$$\psi ' = \dfrac{\phi ' + 2 \frac{\phi}{\eta}}{\eta-\phi}\psi.$$Note that the equations are in normal form. The initial conditions are: $$\psi (1) = \dfrac{\gamma + 1}{\gamma -1},\quad f(1)=\dfrac{2\gamma}{\gamma + 1},\quad \phi (1)=\dfrac{2}{\gamma + 1}.$$The result I'm expecting is more or less this, where density=$\psi$, velocity=$\phi$, pressure=$f$:As you can see, $\phi$ becomes soon linear, with a slope that can be seen to be $\frac{1}{\gamma}$, while $\psi$ and $f$ respectively approach the costant values $0$ and (according to Taylor) $0.436$.  The integration scheme I'm using right now is the midpoint method: $$y_{n+1}=y_n + hF(x_n+\frac{h}{2},y_n+\frac{h}{2}F(x_n,y_n)).$$Everything seems to work quite well until $\eta \approx 0.5$, $\phi'$ seems to approach $0.714=1/\gamma$ (however $f$ at that point is about $0.427$, which is a little lower than its limit and this should not occur). About at that point, $\phi'$ starts again to rise and for $\eta \to 0$ $\phi',\phi,\psi'$ and $\psi$ reach astronomical values. This is a plot of $\phi$ vs. $\eta$; here the ihtegration was done with a step $h=0.0001$:I think that the main problem lies in the approximation errors due to the little numbers that enter the ratios I have to calculate in order to compute $f',\phi ',\psi '$ (that are computed as in the equations above in my program). If it can be of any relevance, I'm writing the program in C.Could you suggest any improvement in the algorithm (or maybe a more robust one) and/or the possible problems which may arise in this computation?  I have tried the substitution $\xi = \phi - \eta/\gamma$ as suggested by Kirill and effectively the system seems to be more stable. However I'm still getting an explosion near $0$: 
Numerical integration algorithm for this set of ODEs
algorithms;numerical analysis
It is entirely possible, although I cannot tell by looking at it, that your system becomes stiff near $\eta=0$. I say this because Mathematica's NDSolve integrates this up to $\eta=3.4\times 10^{-9}$, then complains about a singularity or stiffness. The midpoint method is not A-stable, so that may be the source of inaccuracy: it is unsuitable for stiff problems. The ODE has some kind of singularity near $0$, so maybe the solution blows up because it hits the singularity due to inaccuracy before the method's instability becomes apparent. What do you get if you use a larger time step, like $0.01$ or $0.001$?Try an A-stable method, like the trapezoid rule:$$ y_{n+1} = y_n + \frac{h}{2}(F(x_n,y_n) + F(x_{n+1}, y_{n+1})). $$You'd need to solve the implicit system of nonlinear equations here.Another approach, since you know the asymptotic behaviour of the ODE solution, is to write the unknown functions as perturbations around the known asymptotic behaviour (e.g., $\frac1\gamma\eta+\eta^2\hat\phi(\eta))$, etc., as in a Taylor series), which might make the resulting ODE in $\hat\phi$ easier to solve by (maybe) avoiding the numerical errors in the denominator.