(Sets of) ordinary differential equations (ODE) can be numerically solved using \pstODEsolve from the pst-ode Plain-TeX / LaTeX package.
Integration of the ODEs is done using the Runge-Kutta-Fehlberg method (RKF45).
Run 2 times lualatex or the sequence of commands
latex myfile
dvips myfile
ps2pdf -dNOSAFER myfile.ps
on the TeX input. Either method writes the result table into a text file on the first run. The first method makes use of luapstricks, a PostScript interpreter recently written by Marcel Krüger.
\documentclass{article}
\usepackage{pst-ode}
\usepackage{amsmath}
\usepackage{pgfplots}
\pgfplotsset{compat=newest}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%solve dy/dx=x^2 + y^2 - 1 numerically for different initial values of y in the
%interval x=[-1.1,1.1]; write resulting curves as tables with 100 output points
%into text files
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%y0=-0.5 --> y0=-0.5.dat
\pstODEsolve[algebraicOutputFormat,algebraic,saveData]{%
y0=-0.5% %name of the output file y0=-0.5.dat' }{ t | x[0] %table format iny0=-0.5.dat': x y
}{
-1.1 %integration domain x_min=-1.1
}{
1.1 %integration domain x_max=1.1
}{
100 %number of output points
}{
-0.5 %initial value y0(x_min)
}{
t^2+x[0]^2-1 % right hand side of ODE, note the special notation: x --> t, y --> x[0]
}
%y0=0.0 --> y0=0.0.dat
\pstODEsolve[algebraicOutputFormat,algebraic,saveData]{y0=0.0}{t | x[0]}{-1.1}{1.1}{100}{0.0}{t^2+x[0]^2-1}
%y0=0.5 --> y0=0.5.dat
\pstODEsolve[algebraicOutputFormat,algebraic,saveData]{y0=0.5}{t | x[0]}{-1.1}{1.1}{100}{0.5}{t^2+x[0]^2-1}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
\IfFileExists{y0=-0.5.dat}{}{dummy text\end{document}}
\begin{tikzpicture}
\begin{axis}[
axis equal image, % Unit vectors for both axes have the same length
xmin=-1.1, xmax=1.1, % Axis limits
ymin=-1.1, ymax=1.1,
xtick={-1,-0.5,0,0.5,1}, ytick={-1,-0.5,0,0.5,1}, % Tick marks
no markers,
title={$\dfrac{\mathrm{d}y}{\mathrm{d}x}=x^2+y^2-1$},
view={0}{90},samples=21,domain=-1.1:1.1, y domain=-1.1:1.1, %for direction field
]
\addplot3 [gray, quiver={u={1}, v={x^2+y^2-1}, scale arrows=0.075, every arrow/.append style={-latex}}] (x,y,0);
\addplot table {y0=-0.5.dat};
\addplot table {y0=0.0.dat};
\addplot table {y0=0.5.dat};
\end{axis}
\end{tikzpicture}
\end{document}
