summaryrefslogtreecommitdiff
path: root/macros/latex/contrib/mfpic4ode
diff options
context:
space:
mode:
authorNorbert Preining <norbert@preining.info>2019-09-02 13:46:59 +0900
committerNorbert Preining <norbert@preining.info>2019-09-02 13:46:59 +0900
commite0c6872cf40896c7be36b11dcc744620f10adf1d (patch)
tree60335e10d2f4354b0674ec22d7b53f0f8abee672 /macros/latex/contrib/mfpic4ode
Initial commit
Diffstat (limited to 'macros/latex/contrib/mfpic4ode')
-rw-r--r--macros/latex/contrib/mfpic4ode/README16
-rw-r--r--macros/latex/contrib/mfpic4ode/demo/demo-plain.pdfbin0 -> 213804 bytes
-rw-r--r--macros/latex/contrib/mfpic4ode/demo/demo-plain.tex368
-rw-r--r--macros/latex/contrib/mfpic4ode/demo/demo.pdfbin0 -> 221236 bytes
-rw-r--r--macros/latex/contrib/mfpic4ode/demo/demo.tex381
-rw-r--r--macros/latex/contrib/mfpic4ode/mfpic4ode.dtx429
-rw-r--r--macros/latex/contrib/mfpic4ode/mfpic4ode.ins59
-rw-r--r--macros/latex/contrib/mfpic4ode/mfpic4ode.pdfbin0 -> 231163 bytes
8 files changed, 1253 insertions, 0 deletions
diff --git a/macros/latex/contrib/mfpic4ode/README b/macros/latex/contrib/mfpic4ode/README
new file mode 100644
index 0000000000..12f93d2aa0
--- /dev/null
+++ b/macros/latex/contrib/mfpic4ode/README
@@ -0,0 +1,16 @@
+The mfpic4ode is a set of macros for drawing direction fields, phase
+portraits and trajectories of differential equations and two
+dimensional autonomous systems. The Euler, Runge-Kutta and 4th order
+Runge-Kutta algorithms are available to solve the ODE's. The picture
+is translated into mfpic macros and MetaPost is used to create the
+final drawing. The package is intended for LaTeX, but it can be used
+in plainTeX as well.
+
+Online demonstration of the mfpic4ode macros is available on
+http://wood.mendelu.cz/math/mfpicpreviewer/ as Example 6.
+
+You can get the documentation but compiling mfpic4ode.dtx file and the
+files for installation by compiling mfpic4ode.ins by LaTeX.
+
+The author is supported by the grant FRVS 99/2008 (Fund for
+Developement of Czech Universities). \ No newline at end of file
diff --git a/macros/latex/contrib/mfpic4ode/demo/demo-plain.pdf b/macros/latex/contrib/mfpic4ode/demo/demo-plain.pdf
new file mode 100644
index 0000000000..0d96e43985
--- /dev/null
+++ b/macros/latex/contrib/mfpic4ode/demo/demo-plain.pdf
Binary files differ
diff --git a/macros/latex/contrib/mfpic4ode/demo/demo-plain.tex b/macros/latex/contrib/mfpic4ode/demo/demo-plain.tex
new file mode 100644
index 0000000000..f8e758f3d6
--- /dev/null
+++ b/macros/latex/contrib/mfpic4ode/demo/demo-plain.tex
@@ -0,0 +1,368 @@
+% compile with pdftex demo-plain; mpost portret; pdftex demo-plain
+\input multido
+\input mfpic
+
+\usemetapost
+\opengraphsfile{portret}
+\tlabelsep{3pt}
+\usemplabels
+
+
+\def\frac#1#2{{#1 \over #2}}
+\mfpicunit=1cm
+
+\input mfpic4ode
+
+
+
+
+
+\centerline{\bf Test file for mfpic4ode package}
+\centerline{Robert Ma\v r\'\i k}
+\centerline{January 3, 2008}
+
+\bigskip
+
+See the source file {\tt demo-plain.tex} for comments in the \TeX{} code.
+
+\clipmfpic
+\bigskip
+\centerline{\bf Logistic equation}
+Here we draw a simple picture which describes stability of stationary
+points of teh equation and then draw phase portrait of the equation.
+
+$$ x'=
+ {{r}\cdot\left(1-{x\over K}\right)x}
+$$
+
+% We set parameters for logistic equation
+\mfsrc{r:=1;K:=0.98;}
+
+% We set parameters for drawing and for the ODE solver
+\mfsrc{ODEarrowlength:=0.07; ODEstep:=0.02; ODEstepcount:=500;}
+
+% We define the equation
+\ODEdefineequation{r*y*(1-(y/K))}
+
+Stability and sign of the right--hand side.
+
+\mfpic[5][3]{-0.1}{1.5}{-0.1}{0.5}
+ \axes
+ \xmarks{K}
+ \tlabel[tc](K,0){$K$}
+ \tlabel[bc](0,ypos){$\smash{f(x)}$}
+ \tlabel[cl](xpos,0){$x$}
+
+ % This code draws arrows on x axis, the arrow points to the right if
+ % $f(x)$ is positive and to the left if f(x) is negative. If the
+ % starting point of the arrow is x_0 and the final point x_1 and if
+ % the function changes sign between x_0 and x_1, the arrow from x_0
+ % is not drawn
+ \multido{\r=0.01+0.2}{8}{\ODEharrow{\r}}
+
+ \pen{1pt}
+ \draw[rgb(0,0.5,0)]\parafcn{0,1.5,0.1}{(t,ODErhs(1,t))}
+ \draw[red]\parafcn{0,1.5,0.1}{(t,ODErhs(1,t))}
+ \draw[blue]\parafcn{0,K,0.1}{(t,ODErhs(1,t))}
+\endmfpic
+
+
+\mfsrc{ODEarrowlength:=0.3;}
+
+
+Phase portrait
+
+\mfpic[1][4]{-0.1}{10}{-0.1}{1.5}
+ \axes
+ \ymarks{K}
+ \tlabel[cr](0,K){$K$}
+ \penwd{1pt}
+ \tlabel[tc](xpos,0){$t$}
+ \tlabel[bc](0,ypos){$x$}
+ \pen{0.7pt}
+
+ % This code draws arrows in 19x15 points and draws three integral
+ % curves
+
+ \multido{\r=0.0+0.1}{15}
+ {\multido{\R=0.0+0.5}{19}
+ {\ODEarrow{\R}{\r}}}
+ \trajectories{0,0.3;0,0.01;0,1.4}
+
+\endmfpic
+
+\bigskip
+\centerline{\bf Logistic equation with harvesting}
+
+Similar to the previous picture, but both pictures are drawn together
+to see the relations between them.
+
+$$ x'=
+ {{r}\cdot\left(1-\frac x{{K}}\right)x}-{p}
+$$
+
+% We set parameters for logistic equation with harvesting and define
+% this equation.
+\mfsrc{r:=1;K:=0.98;lov:=0.15;}
+\ODEdefineequation{r*y*(1-(y/K))-lov}
+
+% If the equation possesses stationary points, we store them into
+% variables meza and mezb. If not, we set these variables to negative
+% values (and the are not drawn in view of mfpicclip option.
+\mfsrc{if (r**2*(K**2)-4*r*lov*K)<0: meza:=-1;mezb:=-1.1
+ else: meza:=(r*K-sqrt(r**2*(K**2)-4*r*lov*K))/(2*r);
+ mezb:=(r*K+sqrt(r**2*(K**2)-4*r*lov*K))/(2*r)
+ fi;}
+
+
+\hbox to \hsize{\hss\mfpic[1][4]{-0.1}{10}{-0.1}{1.3}
+ \axes
+ \ymarks{K}
+ \tlabel[cr](0,K){$K$}
+ \penwd{1pt}
+ \tlabel[tc](xpos,0){$t$}
+ \tlabel[bc](0,ypos){$x$}
+
+ % here we draw lines - stationary solutions stored in meza and mezb
+ % variables
+ \draw[gray(0.5)]\dashed\lines{(0,meza),(xpos,meza)}
+ \draw[gray(0.5)]\dashed\lines{(0,mezb),(xpos,mezb)}
+
+ % We draw direction field using metapost cycle. Another option is to
+ % use multido command as in the previous example.
+ \mfsrc{for j=0 step 0.07 until 1.2:
+ for i:=0 step 0.5 until 10:}
+ \ODEarrow{i}{j}
+ \mfsrc{endfor;endfor;}
+
+ % We draw trajectories using black color.
+ \drawcolor{black}
+ \trajectories{0,0.1;0,0.2;0,0.4;0,0.6;0,0.8;0,1.1}
+\endmfpic\qquad
+% On the right we draw the right hand side of the equation
+\mfpic[3][4]{-0.15}{0.6}{-0.1}{1.3}
+ \axes
+ \ymarks{K}
+ \tlabel[cr](0,K){$K$}
+ \tlabel[br](xpos,0){$f(x)$}
+ \tlabel[bc](0,ypos){$x$}
+ \pen{1pt}
+
+ % This code draws the graph of right-hand side of logistic equation
+ % without harvesting.
+ \drawcolor{gray(0.7)}\draw\parafcn{0,1.2,0.1}{(ODErhs(1,t)+lov,t)}
+
+ % This code draws the graph of right-hand side. We use the blue
+ % color for positivce and red color for negative parts. We draw also
+ % arrows which are determined by the sigh of the right-hand side of
+ % the equation.
+ \drawcolor{red}\draw\parafcn{0,1.2,0.1}{(ODErhs(1,t),t)}
+ \drawcolor{blue}\parafcn{meza,mezb,0.05}{(ODErhs(1,t),t)}
+ \multido{\r=0.1+0.15}{7}{\ODEvarrow{\r}}
+\endmfpic\hss}
+
+\bigskip
+\centerline{\bf Three numerical methods for ODEs}
+
+Here we draw solution of ODE using all three available numerical
+methods. We use big step to see the difference between Euler,
+Runge--Kutta and fourth order Runge--Kutta method.
+$$ y'=x+y^3\qquad y(0)=1
+$$
+ \leavevmode
+ \mfpic[20][5]{0}{0.5}{0.9}{2.4}
+ % We set-up parameters
+ \mfsrc{ODEarrowlength:=0.5;}
+ \ODEdefineequation{x+(y**3)}
+ \pen{1pt}\tlabelsep{1pt}
+
+ % We set up parameters for small step
+ \nomplabels
+ \drawcolor{green}
+ \mfsrc{ODEstep:=0.02; ODEstepcount:=30;}
+ \trajectoryRKF{0}{1}
+ \tlabel[tr](0.39,2.4){\bf{Exact solution}}
+
+ % We use bigger step to see the difference between various
+ % methods.
+ \mfsrc{ODEstep:=0.2; ODEstepcount:=2;}
+
+ % We draw trajectory by Euler method
+ \drawcolor{black}
+ \trajectory{0}{1}
+ \tlabel[bl](0.4,1.6){\bf{Euler}}
+ \tlabel[tl](0,1){\bf{$k_1$}}
+ \tlabel[tl](0.2,1.2){\bf{$k_1$ for second step}}
+
+ % We draw trajectory by Runge-Kutta method
+ \drawcolor{rgb(0.5,0.5,0.5)}
+ \trajectoryRK{0}{1}
+ \tlabel[cl](0.4,2.05){\bf{RK}}
+ \tlabel[tl](0.1,1.1){\bf{$k_2$}}
+
+ % We draw trajectory by fourth order Runge-Kutta method
+ \drawcolor{rgb(1,0,0)}
+ \trajectoryRKF{0}{1}
+ \tlabel[tl](0.4,2.15){\bf{RK4}}
+
+ \tlabelsep{3pt}
+
+ % We draw direction field using blue arrows and metapost cycle
+ \penwd{1pt}
+ \drawcolor{blue}\headcolor{blue}
+ \mfsrc{for j=0.9 step 0.1 until 2.3:
+ for i:=0 step 0.05 until 0.5:}
+ \ODEarrow{i}{j}
+ \mfsrc{endfor;endfor;}
+
+
+ \drawcolor{black}
+ \doaxes{lbrt}
+ \bmarks{0,0.2,0.4}
+ \tmarks{0,0.2,0.4}
+ \lmarks{1,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,2,2.2}
+ \pointcolor{red}
+ \pointfilltrue
+ \point[4pt]{(0,1)}
+
+ % We draw some other and bigger arrows in the direction field. The
+ % slopes of hese arrows are important for the first step by Euler
+ % and Ruge-Kutta method and the second step by Euler method.
+ \pen{2pt}
+ \mfsrc{ODEarrowlength:=1;}
+ \colorODEarrowfalse
+ \drawcolor{red}\headcolor{red}
+ \ODEarrows{0,1;0.1,1.1}
+ \ODEarrow{0.2}{1.2}
+
+ \axislabels{b}[tc]{{$0$}0,{$0.2$}0.2,{$0.4$}0.4}
+ \axislabels{l}[cr]{{$0.8$}0.8,{$1$}1,{$1.2$}1.2,{$1.4$}1.4,{$1.6$}1.6,{$1.8$}1.8,{$2$}2,{$2.2$}2.2,{$2.4$}2.4}
+ \endmfpic
+
+
+\bigskip
+\centerline{\bf Autonomous systems}
+
+% The color arrows have no sense in the phase portrait of autonomous
+% system.
+\colorODEarrowfalse
+
+We draw the phase portrait of autonomous system, nulclines, invariant
+set between nulclines, trajectories. We draw arrows in regular grid
+and add few more arrows on nulclines and outside the regular grid.
+
+\mfsrc{TIMEstep:=0.05; TIMEend:=30;}
+
+\mfpic
+ [0.5]{-2}{15}{-2}{15}
+ \nomplabels
+ \tlabel[cc](8,15.5){\bf Competing species}
+ \usemplabels
+
+ % We set up parameters, define equations, define functions which
+ % describe nulclines and store stationary points into variables z10,
+ % z11, z12, z1.
+
+ \mfsrc{a:=11;b:=1;c:=0.8;k:=10;l:=1.1;m:=1.2;}
+ \mfsrc{ODEarrowlength:=0.3;}
+ \ASdefineequations{x*(a-b*x-c*y)}{y*(k-l*x-m*y)}
+ \fdef{xnulklina}{x}{(a-b*x)/c}
+ \fdef{ynulklina}{x}{(k-l*x)/m}
+ \mfsrc{z10=(0,a/c);z11=(0,k/m);z12=(a/b,0);z13=(k/l,0);}
+
+ % Here we draw a gray polygon - invariant set fot the system.
+ \pen{0.3pt}
+ \gfill[gray(0.7)]\lclosed\lines{z10,z11,z13,z12}
+ \axes
+ \tlabel[bc](0,ypos){$y$}
+ \tlabel[cl](xpos,0){$x$}
+
+ \pointsize=3pt
+
+ \pointfilltrue\pointcolor{red}\point{(a/b,0)}
+ \draw[red]\function{0,a/b,1}{xnulklina(x)}
+ \draw[red]\lines{(0,0),(0,ypos)}
+ \tlabel[cr](x10,y10){$a\over c$}
+ \tlabel[tc](x12,y12){$a\over b$}
+
+ \pointfilltrue\pointcolor{green}\point{(0,k/m)}
+ \draw[green]\function{0,k/l,1}{ynulklina(x)}
+ \draw[green]\lines{(0,0),(xpos,0)}
+ \tlabel[cr](x11,y11){$\alpha\over \beta$}
+ \tlabel[tc](x13,y13){$\alpha\over \gamma$}
+
+
+ \penwd{1.5pt}
+ \drawcolor{gray(0.25)}\headcolor{gray(0.25)}
+ \ASarrows{0,6;0,11;7,0;9,0;4.5,0;13,0;2,0;0,14;0,2}
+ \ASarrows{3,3;8,7;13,2;3,14}
+ \ASarrows{4,ynulklina(4);5,ynulklina(5);6,ynulklina(6);
+ 1.7,ynulklina(1.7)}
+ \ASarrows{0.5,xnulklina(0.5);1.8,xnulklina(1.8);6,xnulklina(6);4,xnulklina(4);
+ 7.5,xnulklina(7.5)}
+
+ \drawcolor{black}
+ \AStrajectories{12,12;8,12;4,12;0.2,3;1,3;6,12;1,1;1,0.12;0.3,14}
+ \mfsrc{TIMEstep:=-0.05; TIMEend:=5;}
+ \AStrajectories{12,12;8,12;4,12;0.2,3;1,3;6,12;1,1;1,0.12;0.3,14}
+ \penwd{0.5pt}
+ \drawcolor{gray(0.5)}
+ \headcolor{gray(0.5)}
+ \multido{\r=0.5+1}{15}{\multido{\R=0.5+1}{15}{\ASarrow{\R}{\r}}}
+
+
+\endmfpic
+
+
+\break
+\bigskip
+\centerline{\bf Predator prey system with HollingII response function}
+
+
+
+
+\mfpic[2]{-0.1}{4}{-0.1}{3}
+ % we define functions and parameters, right hand sides of the system
+ % and a function which defines nulcline
+ \mfsrc{r:=1;K:=3;a:=1;k:=0.8;P:=1;Alfa:=0.42;}
+ \fdef{funkceV}{x}{a*x/(x+P)}
+ \ASdefineequations{r*x*(1-(x/K))-funkceV(x)*y}{(-Alfa+k*funkceV(x))*y}
+ \fdef{xnulklina}{x}{r*(1-(x/K))*(x+P)/a}
+ \mfsrc{ODEarrowlength:=0.2;}
+
+ % Here we draw axes and nulclines
+ \axes
+ \tlabel[bc](0,ypos){$y$}
+ \tlabel[cl](xpos,0){$x$}
+ \draw[red]\function{0,K,0.1}{xnulklina(x)}
+ \draw[green]\lines{(P/((k*a/Alfa)-1),0),(P/((k*a/Alfa)-1),ypos)}
+
+ % Here we draw some arrows on nulclines and then arrows in the plane
+ \penwd{0.5pt}
+ \drawcolor{gray(0.25)}\headcolor{gray(0.25)}
+ \ASarrows{P/((k*a/Alfa)-1),1;P/((k*a/Alfa)-1),2;P/((k*a/Alfa)-1),0.5;P/((k*a/Alfa)-1),1.5}
+ \ASarrows{0,xnulklina(0);0.5,xnulklina(0.5);1,xnulklina(1);1.5,xnulklina(1.5);2,xnulklina(2);2.25,xnulklina(2.25)}
+ \multido{\r=0.1+0.25}{20}{\multido{\R=0.1+0.25}{20}{\ASarrow{\R}{\r}}}
+
+ %% We draw trajectory with IC x=2, y=2
+ \drawcolor{black}
+ \mfsrc{TIMEstep:=0.05; TIMEend:=10;}
+ \AStrajectoryRKF{2}{2}
+ %% We continue the trajectory (spiral) from the last point
+ \AStrajectoryRKF{x1}{y1}
+ \AStrajectoryRKF{x1}{y1}
+ \AStrajectoryRKF{x1}{y1}
+ \AStrajectoryRKF{x1}{y1}
+ \AStrajectoryRKF{x1}{y1}
+ \AStrajectoryRKF{x1}{y1}
+ \AStrajectoryRKF{x1}{y1}
+
+ %% We continue backwards
+ \mfsrc{TIMEstep:=-0.05;}
+ \AStrajectoryRKF{2}{2}
+
+\endmfpic
+
+
+\closegraphsfile
+\end
diff --git a/macros/latex/contrib/mfpic4ode/demo/demo.pdf b/macros/latex/contrib/mfpic4ode/demo/demo.pdf
new file mode 100644
index 0000000000..256f443f0b
--- /dev/null
+++ b/macros/latex/contrib/mfpic4ode/demo/demo.pdf
Binary files differ
diff --git a/macros/latex/contrib/mfpic4ode/demo/demo.tex b/macros/latex/contrib/mfpic4ode/demo/demo.tex
new file mode 100644
index 0000000000..0a56345461
--- /dev/null
+++ b/macros/latex/contrib/mfpic4ode/demo/demo.tex
@@ -0,0 +1,381 @@
+% compile with pdflatex demo; mpost portret; pdflatex demo
+\documentclass{article}
+
+\usepackage{multido}
+\usepackage{amsmath}
+%\usepackage{color}
+\usepackage%[clip]
+{mfpic}
+\usemetapost
+\opengraphsfile{portret}
+\tlabelsep{3pt}
+\usemplabels
+
+
+\mfpicunit=1cm
+
+
+
+\usepackage{mfpic4ode}
+
+
+\begin{document}
+
+\title{Test file for mfpic4ode package}
+\author{Robert Ma\v r\'\i k}
+\date\today
+\maketitle
+
+See the source file \texttt{demo.tex} for comments in the \TeX{} code.
+
+\clipmfpic
+\section{Logistic equation}
+Here we draw a simple picture which describes stability of stationary
+points of teh equation and then draw phase portrait of the equation.
+
+\begin{center}
+ \begin{equation}
+ x'=
+ {{r}\cdot\left(1-\frac x{{K}}\right)x}
+ \end{equation}
+\end{center}
+
+% We set parameters for logistic equation
+\mfsrc{r:=1;K:=0.98;}
+
+% We set parameters for drawing and for the ODE solver
+\mfsrc{ODEarrowlength:=0.07; ODEstep:=0.02; ODEstepcount:=500;}
+
+% We define the equation
+\ODEdefineequation{r*y*(1-(y/K))}
+
+Stability and sign of the right--hand side.
+
+\begin{mfpic}[5][3]{-0.1}{1.5}{-0.1}{0.5}
+ \axes
+ \xmarks{K}
+ \tlabel[tc](K,0){$K$}
+ \tlabel[bc](0,ypos){$\smash{f(x)}$}
+ \tlabel[cl](xpos,0){$x$}
+
+ % This code draws arrows on x axis, the arrow points to the right if
+ % $f(x)$ is positive and to the left if f(x) is negative. If the
+ % starting point of the arrow is x_0 and the final point x_1 and if
+ % the function changes sign between x_0 and x_1, the arrow from x_0
+ % is not drawn
+ \multido{\r=0.01+0.2}{8}{\ODEharrow{\r}}
+
+ \pen{1pt}
+ \draw[rgb(0,0.5,0)]\parafcn{0,1.5,0.1}{(t,ODErhs(1,t))}
+ \draw[red]\parafcn{0,1.5,0.1}{(t,ODErhs(1,t))}
+ \draw[blue]\parafcn{0,K,0.1}{(t,ODErhs(1,t))}
+\end{mfpic}
+
+
+\mfsrc{ODEarrowlength:=0.3;}
+
+
+Phase portrait
+
+\begin{mfpic}[1][4]{-0.1}{10}{-0.1}{1.5}
+ \axes
+ \ymarks{K}
+ \tlabel[cr](0,K){$K$}
+ \penwd{1pt}
+ \tlabel[tc](xpos,0){$t$}
+ \tlabel[bc](0,ypos){$x$}
+ \pen{0.7pt}
+
+ % This code draws arrows in 19x15 points and draws three integral
+ % curves
+
+ \multido{\r=0.0+0.1}{15}
+ {\multido{\R=0.0+0.5}{19}
+ {\ODEarrow{\R}{\r}}}
+ \trajectories{0,0.3;0,0.01;0,1.4}
+
+\end{mfpic}
+
+
+\section{Logistic equation with harvesting}
+
+Similar to the previous picture, but both pictures are drawn together
+to see the relations between them.
+
+\begin{center}
+ \begin{equation}
+ x'=
+ {{r}\cdot\left(1-\frac x{{K}}\right)x}-{p}
+ \end{equation}
+\end{center}
+
+% We set parameters for logistic equation with harvesting and define
+% this equation.
+\mfsrc{r:=1;K:=0.98;lov:=0.15;}
+\ODEdefineequation{r*y*(1-(y/K))-lov}
+
+% If the equation possesses stationary points, we store them into
+% variables meza and mezb. If not, we set these variables to negative
+% values (and the are not drawn in view of mfpicclip option.
+\mfsrc{if (r**2*(K**2)-4*r*lov*K)<0: meza:=-1;mezb:=-1.1
+ else: meza:=(r*K-sqrt(r**2*(K**2)-4*r*lov*K))/(2*r);
+ mezb:=(r*K+sqrt(r**2*(K**2)-4*r*lov*K))/(2*r)
+ fi;}
+
+
+\hbox to \hsize{\hss\begin{mfpic}[1][4]{-0.1}{10}{-0.1}{1.3}
+ \axes
+ \ymarks{K}
+ \tlabel[cr](0,K){$K$}
+ \penwd{1pt}
+ \tlabel[tc](xpos,0){$t$}
+ \tlabel[bc](0,ypos){$x$}
+
+ % here we draw lines - stationary solutions stored in meza and mezb
+ % variables
+ \draw[gray(0.5)]\dashed\lines{(0,meza),(xpos,meza)}
+ \draw[gray(0.5)]\dashed\lines{(0,mezb),(xpos,mezb)}
+
+ % We draw direction field using metapost cycle. Another option is to
+ % use multido command as in the previous example.
+ \mfsrc{for j=0 step 0.07 until 1.2:
+ for i:=0 step 0.5 until 10:}
+ \ODEarrow{i}{j}
+ \mfsrc{endfor;endfor;}
+
+ % We draw trajectories using black color.
+ \drawcolor{black}
+ \trajectories{0,0.1;0,0.2;0,0.4;0,0.6;0,0.8;0,1.1}
+\end{mfpic}\qquad
+% On the right we draw the right hand side of the equation
+\begin{mfpic}[3][4]{-0.15}{0.6}{-0.1}{1.3}
+ \axes
+ \ymarks{K}
+ \tlabel[cr](0,K){$K$}
+ \tlabel[br](xpos,0){$f(x)$}
+ \tlabel[bc](0,ypos){$x$}
+ \pen{1pt}
+
+ % This code draws the graph of right-hand side of logistic equation
+ % without harvesting.
+ \drawcolor{gray(0.7)}\draw\parafcn{0,1.2,0.1}{(ODErhs(1,t)+lov,t)}
+
+ % This code draws the graph of right-hand side. We use the blue
+ % color for positivce and red color for negative parts. We draw also
+ % arrows which are determined by the sigh of the right-hand side of
+ % the equation.
+ \drawcolor{red}\draw\parafcn{0,1.2,0.1}{(ODErhs(1,t),t)}
+ \drawcolor{blue}\parafcn{meza,mezb,0.05}{(ODErhs(1,t),t)}
+ \multido{\r=0.1+0.15}{7}{\ODEvarrow{\r}}
+\end{mfpic}\hss}
+
+\section{Three numerical methods for ODEs}
+
+Here we draw solution of ODE using all three available numerical
+methods. We use big step to see the difference between Euler,
+Runge--Kutta and fourth order Runge--Kutta method.
+\begin{center}
+ \begin{equation*}
+ \boxed{\begin{aligned}y'&=x+y^3\\y(0)&=1
+ \end{aligned}}\qquad
+ \boxed{\begin{aligned}x_{n+1}&=x_n+h\\
+ y_{n+1}&=y_n+kh\\h&=0.2
+ \end{aligned}}
+ \end{equation*}
+ \leavevmode
+ \begin{mfpic}[20][5]{0}{0.5}{0.9}{2.4}
+ % We set-up parameters
+ \mfsrc{ODEarrowlength:=0.5;}
+ \ODEdefineequation{x+(y**3)}
+ \pen{1pt}\tlabelsep{1pt}
+
+ % We set up parameters for small step
+ \nomplabels
+ \drawcolor{green}
+ \mfsrc{ODEstep:=0.02; ODEstepcount:=30;}
+ \dashed\trajectoryRKF{0}{1}
+ \tlabel[tr](0.39,2.4){\bf{Exact solution}}
+
+ % We use bigger step to see the difference between various
+ % methods.
+ \mfsrc{ODEstep:=0.2; ODEstepcount:=2;}
+
+ % We draw trajectory by Euler method
+ \drawcolor{black}
+ \trajectory{0}{1}
+ \tlabel[bl](0.4,1.6){\bf{Euler}}
+ \tlabel[tl](0,1){\bf{$k_1$}}
+ \tlabel[tl](0.2,1.2){\bf{$k_1$ for second step}}
+
+ % We draw trajectory by Runge-Kutta method
+ \drawcolor{rgb(0.5,0.5,0.5)}
+ \trajectoryRK{0}{1}
+ \tlabel[cl](0.4,2.05){\bf{RK}}
+ \tlabel[tl](0.1,1.1){\bf{$k_2$}}
+
+ % We draw trajectory by fourth order Runge-Kutta method
+ \drawcolor{rgb(1,0,0)}
+ \trajectoryRKF{0}{1}
+ \tlabel[tl](0.4,2.15){\bf{RK4}}
+
+ \tlabelsep{3pt}
+
+ % We draw direction field using blue arrows and metapost cycle
+ \penwd{1pt}
+ \drawcolor{blue}\headcolor{blue}
+ \mfsrc{for j=0.9 step 0.1 until 2.3:
+ for i:=0 step 0.05 until 0.5:}
+ \ODEarrow{i}{j}
+ \mfsrc{endfor;endfor;}
+
+
+ \drawcolor{black}
+ \doaxes{lbrt}
+ \bmarks{0,0.2,0.4}
+ \tmarks{0,0.2,0.4}
+ \lmarks{1,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,2,2.2}
+ \pointcolor{red}
+ \pointfilltrue
+ \point[4pt]{(0,1)}
+
+ % We draw some other and bigger arrows in the direction field. The
+ % slopes of hese arrows are important for the first step by Euler
+ % and Ruge-Kutta method and the second step by Euler method.
+ \pen{2pt}
+ \mfsrc{ODEarrowlength:=1;}
+ \colorODEarrowfalse
+ \drawcolor{red}\headcolor{red}
+ \ODEarrows{0,1;0.1,1.1}
+ \ODEarrow{0.2}{1.2}
+
+ \axislabels{b}[tc]{{$0$}0,{$0.2$}0.2,{$0.4$}0.4}
+ \axislabels{l}[cr]{{$0.8$}0.8,{$1$}1,{$1.2$}1.2,{$1.4$}1.4,{$1.6$}1.6,{$1.8$}1.8,{$2$}2,{$2.2$}2.2,{$2.4$}2.4}
+ \end{mfpic}
+
+\end{center}
+
+
+\section{Autonomous systems}
+
+% The color arrows have no sense in the phase portrait of autonomous
+% system.
+\colorODEarrowfalse
+
+We draw the phase portrait of autonomous system, nulclines, invariant
+set between nulclines, trajectories. We draw arrows in regular grid
+and add few more arrows on nulclines and outside the regular grid.
+
+\mfsrc{TIMEstep:=0.05; TIMEend:=30;}
+
+\begin{mfpic}
+ [0.5]{-2}{15}{-2}{15}
+ \nomplabels
+ \tlabel[cc](8,15.5){\bf Competing species}
+ \usemplabels
+
+ % We set up parameters, define equations, define functions which
+ % describe nulclines and store stationary points into variables z10,
+ % z11, z12, z1.
+
+ \mfsrc{a:=11;b:=1;c:=0.8;k:=10;l:=1.1;m:=1.2;}
+ \mfsrc{ODEarrowlength:=0.3;}
+ \ASdefineequations{x*(a-b*x-c*y)}{y*(k-l*x-m*y)}
+ \fdef{xnulklina}{x}{(a-b*x)/c}
+ \fdef{ynulklina}{x}{(k-l*x)/m}
+ \mfsrc{z10=(0,a/c);z11=(0,k/m);z12=(a/b,0);z13=(k/l,0);}
+
+ % Here we draw a gray polygon - invariant set fot the system.
+ \pen{0.3pt}
+ \gfill[gray(0.7)]\lclosed\lines{z10,z11,z13,z12}
+ \axes
+ \tlabel[bc](0,ypos){$y$}
+ \tlabel[cl](xpos,0){$x$}
+
+ \pointsize=3pt
+
+ \pointfilltrue\pointcolor{red}\point{(a/b,0)}
+ \draw[red]\function{0,a/b,1}{xnulklina(x)}
+ \draw[red]\lines{(0,0),(0,ypos)}
+ \tlabel[cr](x10,y10){$a\over c$}
+ \tlabel[tc](x12,y12){$a\over b$}
+
+ \pointfilltrue\pointcolor{green}\point{(0,k/m)}
+ \draw[green]\function{0,k/l,1}{ynulklina(x)}
+ \draw[green]\lines{(0,0),(xpos,0)}
+ \tlabel[cr](x11,y11){$\alpha\over \beta$}
+ \tlabel[tc](x13,y13){$\alpha\over \gamma$}
+
+
+ \penwd{1.5pt}
+ \drawcolor{gray(0.25)}\headcolor{gray(0.25)}
+ \ASarrows{0,6;0,11;7,0;9,0;4.5,0;13,0;2,0;0,14;0,2}
+ \ASarrows{3,3;8,7;13,2;3,14}
+ \ASarrows{4,ynulklina(4);5,ynulklina(5);6,ynulklina(6);
+ 1.7,ynulklina(1.7)}
+ \ASarrows{0.5,xnulklina(0.5);1.8,xnulklina(1.8);6,xnulklina(6);4,xnulklina(4);
+ 7.5,xnulklina(7.5)}
+
+ \drawcolor{black}
+ \AStrajectories{12,12;8,12;4,12;0.2,3;1,3;6,12;1,1;1,0.12;0.3,14}
+ \mfsrc{TIMEstep:=-0.05; TIMEend:=5;}
+ \AStrajectories{12,12;8,12;4,12;0.2,3;1,3;6,12;1,1;1,0.12;0.3,14}
+ \penwd{0.5pt}
+ \drawcolor{gray(0.5)}
+ \headcolor{gray(0.5)}
+ \multido{\r=0.5+1}{15}{\multido{\R=0.5+1}{15}{\ASarrow{\R}{\r}}}
+
+
+\end{mfpic}
+
+
+\begin{center}
+ \textbf{Pedator prey system with HollingII response function}
+\end{center}
+
+
+
+\begin{mfpic}[2]{-0.1}{4}{-0.1}{3}
+ % we define functions and parameters, right hand sides of the system
+ % and a function which defines nulcline
+ \mfsrc{r:=1;K:=3;a:=1;k:=0.8;P:=1;Alfa:=0.42;}
+ \fdef{funkceV}{x}{a*x/(x+P)}
+ \ASdefineequations{r*x*(1-(x/K))-funkceV(x)*y}{(-Alfa+k*funkceV(x))*y}
+ \fdef{xnulklina}{x}{r*(1-(x/K))*(x+P)/a}
+ \mfsrc{ODEarrowlength:=0.2;}
+
+ % Here we draw axes and nulclines
+ \axes
+ \tlabel[bc](0,ypos){$y$}
+ \tlabel[cl](xpos,0){$x$}
+ \draw[red]\function{0,K,0.1}{xnulklina(x)}
+ \draw[green]\lines{(P/((k*a/Alfa)-1),0),(P/((k*a/Alfa)-1),ypos)}
+
+ % Here we draw some arrows on nulclines and then arrows in the plane
+ \penwd{0.5pt}
+ \drawcolor{gray(0.25)}\headcolor{gray(0.25)}
+ \ASarrows{P/((k*a/Alfa)-1),1;P/((k*a/Alfa)-1),2;P/((k*a/Alfa)-1),0.5;P/((k*a/Alfa)-1),1.5}
+ \ASarrows{0,xnulklina(0);0.5,xnulklina(0.5);1,xnulklina(1);1.5,xnulklina(1.5);2,xnulklina(2);2.25,xnulklina(2.25)}
+ \multido{\r=0.1+0.25}{20}{\multido{\R=0.1+0.25}{20}{\ASarrow{\R}{\r}}}
+
+ %% We draw trajectory with IC x=2, y=2
+ \drawcolor{black}
+ \mfsrc{TIMEstep:=0.05; TIMEend:=10;}
+ \AStrajectoryRKF{2}{2}
+ %% We continue the trajectory (spiral) from the last point
+ \AStrajectoryRKF{x1}{y1}
+ \AStrajectoryRKF{x1}{y1}
+ \AStrajectoryRKF{x1}{y1}
+ \AStrajectoryRKF{x1}{y1}
+ \AStrajectoryRKF{x1}{y1}
+ \AStrajectoryRKF{x1}{y1}
+ \AStrajectoryRKF{x1}{y1}
+
+ %% We continue backwards
+ \mfsrc{TIMEstep:=-0.05;}
+ \AStrajectoryRKF{2}{2}
+
+\end{mfpic}
+
+
+
+\end{document}
diff --git a/macros/latex/contrib/mfpic4ode/mfpic4ode.dtx b/macros/latex/contrib/mfpic4ode/mfpic4ode.dtx
new file mode 100644
index 0000000000..f8e337b3df
--- /dev/null
+++ b/macros/latex/contrib/mfpic4ode/mfpic4ode.dtx
@@ -0,0 +1,429 @@
+% \iffalse meta-comment
+%
+% Copyright (C) 2006-2008 by Robert Marik <marik@mendelu.cz>
+% ----------------------------------------------------------
+%
+% This file may be distributed and/or modified under the
+% conditions of the LaTeX Project Public License, either version 1.2
+% of this license or (at your option) any later version.
+% The latest version of this license is in:
+%
+% http://www.latex-project.org/lppl.txt
+%
+% and version 1.2 or later is part of all distributions of LaTeX
+% version 1999/12/01 or later.
+%
+% \fi
+%
+% \iffalse
+%<*driver>
+\ProvidesFile{mfpic4ode.dtx}
+%</driver>
+%<sty>\NeedsTeXFormat{LaTeX2e}[1999/12/01]
+%<sty>\ProvidesPackage{mfpic4ode}
+%<*sty>
+ [2009/04/15 v0.3 mfpic4ode.dtx file]
+%</sty>
+%
+%<*driver>
+\documentclass{ltxdoc}
+\EnableCrossrefs
+\CodelineIndex
+\RecordChanges
+\begin{document}
+ \DocInput{mfpic4ode.dtx}
+ \PrintChanges
+ \PrintIndex
+\end{document}
+%</driver>
+% \fi
+%
+% \CheckSum{185}
+%
+% \CharacterTable
+% {Upper-case \A\B\C\D\E\F\G\H\I\J\K\L\M\N\O\P\Q\R\S\T\U\V\W\X\Y\Z
+% Lower-case \a\b\c\d\e\f\g\h\i\j\k\l\m\n\o\p\q\r\s\t\u\v\w\x\y\z
+% Digits \0\1\2\3\4\5\6\7\8\9
+% Exclamation \! Double quote \" Hash (number) \#
+% Dollar \$ Percent \% Ampersand \&
+% Acute accent \' Left paren \( Right paren \)
+% Asterisk \* Plus \+ Comma \,
+% Minus \- Point \. Solidus \/
+% Colon \: Semicolon \; Less than \<
+% Equals \= Greater than \> Question mark \?
+% Commercial at \@ Left bracket \[ Backslash \\
+% Right bracket \] Circumflex \^ Underscore \_
+% Grave accent \` Left brace \{ Vertical bar \|
+% Right brace \} Tilde \~}
+%
+%
+% \changes{v0.2}{2008/02/03}{First public version}
+% \changes{v0.3}{2009/04/15}{Added connect environment to paths}
+% \changes{v0.4}{2010/04/07}{Updated documentations}
+%
+% \GetFileInfo{mfpic4ode.dtx}
+%
+% \DoNotIndex{\newcommand,\newenvironment}
+%
+% \def\fileversion{0.4}
+% \def\filedate{2010/04/07}
+%
+% \title{The \textsf{mfpic4ode} package\thanks{This document
+% corresponds to \textsf{mfpic4ode}~\fileversion, dated
+% \filedate.}} \author{Robert Marik \\ \texttt{marik@mendelu.cz}
+% \thanks{Supported by the grant FRV\v{S} 99/2008 (Fund for
+% Developement of Czech Universities).}}
+% \maketitle
+%
+% \section{Introduction}
+%
+% The package |mfpic4ode| is a set of macros for drawing phase
+% portraits and integral curves of differential equations and
+% autonomous systems using |mfpic| macros. These macros have been used
+% by the author to prepare some pictures for classrooms and the
+% results seem to be acceptable for this purpose, but always remember
+% that due to the fixed points arithmetics in Metapost, the error in
+% computations could be significant. Another excellent tool which can
+% be used to draw trajectories is Sage\TeX{} which gives you full
+% power of computer algebra system Sage in \LaTeX.
+%
+% \section{Usage}
+% You can load the package in \LaTeX{} using standard
+% |\usepackage{mfpic4ode}| command, or you can use the macros in
+% plain\TeX{} and load by |\input mfpic4ode.tex| command.
+%
+% \subsection{First order differential equation}
+%
+% To draw phase portrait of first order ordinary differential equation
+% $$
+% y'=f(x,y)
+% $$
+% we define commands |\ODEarrow| for drawing element of direction
+% field and |\trajectory|, |\trajectoryRK| and |\trajectoryRKF| for
+% drawing integral curves using Euler, second order Runge-Kutta and
+% fourth order Runge-Kutta methods, respectively. Some important
+% parameters, such as the number of steps, the length of step or the
+% function from the right-hand side of the equations are stored in
+% MetaPost variables and to keep the package simple and short, these
+% variables are accessible using |\mfsrc| command.
+%
+% \DescribeMacro{\ifcolorODEarrow} \DescribeMacro{\colorODEarrowtrue}
+% \DescribeMacro{\colorODEarrowfalse} If the \TeX{} boolean variable
+% |\ifcolorODEarrow| is true, then the arrows from direction field are
+% blue if the solution is increasing and red if decreasing. If
+% |\ifcolorODEarrow| is false, the mfpic color from |\drawcolor| and
+% |\headcolor| macros is used. More precisely,
+% \begin{itemize}
+% \item if |\ifcolorODEarrow| is true and $f(x_0,y_0)>0$, then the
+% arrow at the point $(x_0,y_0)$ is blue
+% \item if |\ifcolorODEarrow| is true and $f(x_0,y_0)\leq 0$, then the
+% arrow at the point $(x_0,y_0)$ is red
+% \item if |\ifcolorODEarrow| is false, then color from |\drawcolor|
+% is used to draw the body of an arrow and color |\headcolor| is
+% used to draw the head.
+% \end{itemize}
+% Arrows are drawn using mfpic |\draw\arrow\lines{...}| command and
+% hence the parameters for customizing shape and size of the head from
+% mfpic are also available. The MetaPost variable |ODEarrowlength| is
+% used to customize the length of each arrow. If the arrow is
+% horizontal, then the length of the arrow in mfpic coordinates is
+% equal to |ODEarrowlength/xscale|. (This fixes the case when
+% different |xscale| and |yscale| are used. All arrows have the same
+% length.) You can set this variable using |\mfsrc| command, you can
+% write e.g.
+% \begin{verbatim}
+% \mfsrc{ODEarrowlength:=0.07;}
+% \end{verbatim}
+% in your document.
+%
+% To draw arrows in regular rectangular grid you should use the
+% |\ODEarrow| macro in a double cycle such as
+% \begin{verbatim}
+% \mfsrc{for j=0 step 0.07 until 1.2: for i:=0 step 0.5 until 10:}
+% \ODEarrow{i}{j}
+% \mfsrc{endfor;endfor;}
+% \end{verbatim}
+% or using the |multido| package
+% \begin{verbatim}
+% \multido{\r=0.0+0.1}{15}{\multido{\R=0.0+0.5}{19}{\ODEarrow{\R}{\r}}}
+% \end{verbatim}
+%
+% \DescribeMacro{\ODEdefineequation{f(x,y)}} The macro |\ODEdefineequation| is
+% used to save the right hand side of the ODE, i. e. the function
+% $f(x,y)$. You should write the expression in the MetaPost format,
+% the independent variable is supposed to be $x$, the dependent
+% variable is $y$.
+%
+% \DescribeMacro{\trajectory\{x0\}\{y0\}}
+% \DescribeMacro{\trajectoryRK\{x0\}\{y0\}}
+% \DescribeMacro{\trajectoryRKF\{x0\}\{y0\}} The macros |\trajectory|,
+% |\trajectoryRK| and |\trajectoryRKF| are used for drawing integral
+% curves with initial condition $y(x_0)=y_0$ using Euler, second order
+% Runge-Kutta and fourth order Runge-Kutta methods, respectively.
+% The length of each step is stored in MetaPost variable
+% |ODEstepcount|, the length of each step is in the MetaPost variable
+% |ODEstep|. You can set these variables using |\mfsrc| macro as
+% follows
+% \begin{verbatim}
+% \mfsrc{ODEstep:=0.02; ODEstepcount:=500;}
+% \end{verbatim}
+% The integral curve is drawn from short linear parts using |\ODEline|
+% command which expands to |\lines| command from |mfpic| package by
+% default. These linear parts are connected in connect environment and
+% this allows to use prefixes like |\dotted| or |\dashed| to the
+% trajectories. A simple test is used to keep the arithmetics in
+% reasonable bounds: if after the step the curve leaves the horizontal
+% strip between |yneg| and |ypos| variables, then the evaluation is
+% stopped (in fact, in this case we do not change the independent
+% variable and we do the remaining steps with the same last point).
+% Recall that |yneg| and |ypos| variables are set when you call
+% |mfpicture| environment. If you call the environment as follows
+% \begin{verbatim}
+% \begin{mfpic}[5][3]{-0.1}{1.5}{-0.1}{0.5}
+% ...........
+% \end{mfpic}
+% \end{verbatim}
+% then no more than one short linear part of the integral curve is
+% outside the horizontal strip between $y=-0.1$ and $y=0.5$.
+%
+% \DescribeMacro{\trajectories}\DescribeMacro{\ODEarrows} To draw more
+% trajectories you can use |\trajectories| command. The command
+% |\trajectories{x1,y1;x2,y2;x3,y3;....;xn,yn}| expands to $n$
+% |\trajectoryRKF| commands with initial conditions $y(x_i)=y_i$ for
+% $i=1..n$. In a similar way,
+% |\ODEarrows{x1,y1;x2,y2;x3,y3;....;xn,yn}| expands into $n$
+% |\ODEarrow| commands.
+%
+% \subsection{Two-dimensional autonomous systems}
+% Trajectories for two-dimensional autonomous system
+% \begin{eqnarray*}
+% x'&=&f(x,y)\\y'&=&g(x,y)
+% \end{eqnarray*}
+% are drawn using a very simple method based on the direction field.
+% This could be improved in the next release of the package, but till
+% now the results obtained in this way are qualitatively correct and
+% sufficiently accurate (remember that you cannot expect accurate
+% approximation due to the limitation of arithmetics in MetaPost).
+% Anyway, some users may prefer the fourth order Runge--Kutta method.
+%
+% The macros |\ASdefineequation| |\ASarrow|, |\AStrajectory|,
+% |\AStrajectoryRKF|, |\ASarrows| and |\AStrajectories| are
+% counterparts to their |\ODE....| versions. The last point of each
+% trajectory is stored in the |x1| and |x2| MetaPost variables. Hence,
+% you can say |\AStrajectory{2}{2}| to draw trajectory with initial
+% conditions $x(0)=2$, $y(0)=2$ and then you can continue this
+% trajectory using |\AStrajectory{x1}{y1}| command. The macro
+% |\AStrajectory| uses |ODEstep| and |ODEstepcount| variables, the
+% macro |\AStrajectoryRKF| uses |TIMEstep| and |TIMEend| variables do
+% perform the steps in the numerical solution. The number of steps is
+% in the latter case evaluated as absolute value of the quotient
+% |TIMEend/TIMEstep|. You can use negative value for |TIMEstep| to
+% continue the trajectory backwards.
+%
+% \section{Troubleshooting}
+% \subsection{The catcode of @ is messed} We set the category of @ to
+% 11 (letter) when we load the package and at the end of definitions
+% for mfpic4ode we set the category to 12. This could be a source of
+% rare problems, if you use different value in your document.
+%
+% \subsection{Metapost: Not implemented: (unknown numeric) \dots }
+% You have to set |ODEstep|, |ODEstepcount|, |TIMEstep| and |TIMEend|
+% other variables using |\mfsrc| command (depending on the type of the
+% problem).
+%
+% \StopEventually{}
+%
+% \section{Implementation}
+% \begin{macrocode}
+%<*tex>
+\catcode`\@=11
+
+\newif\ifcolorODEarrow
+%%%\colorODEarrowfalse
+\colorODEarrowtrue
+
+%%% The line from one point to another
+\def\ODEline#1#2{\lines{#1,#2}}
+
+%%% The variable ODErhs is used to store the function from the right
+%%% hand side of ODE in the form y'=f(x,y). We use command
+%%% ODEdefineequation to set up this variable.
+\def\ODEdefineequation#1{\fdef{ODErhs}{x,y}{#1}}
+
+%%% Integral curve using Euler method. The step of this method is
+%%% ODEstep and the number of steps is ODEstepcount. The points are
+%%% stored in metapost variables x1,y1.
+\def\trajectory#1#2{
+ \begin{connect}
+ \mfsrc{x1:=#1;y1:=#2;
+ for i=1 upto ODEstepcount:
+ x2:=x1+ODEstep;
+ y2:=y1+ODEstep*ODErhs(x1,y1);}
+ \ODEline{z1}{z2}
+ \mfsrc{
+ if ((y2>yneg) and (y2<ypos)): x1:=x2; y1:=y2 fi;
+ endfor
+ }
+ \end{connect}
+}
+
+%%% Integral curve using Runge--Kutta method.
+\def\trajectoryRK#1#2{
+ \begin{connect}
+ \mfsrc{x1:=#1;y1:=#2;
+ for i=1 upto ODEstepcount:
+ k1:=ODErhs(x1,y1);
+ x3:=x1+(ODEstep/2);
+ y3:=y1+k1*(ODEstep/2);
+ k2:=ODErhs(x3,y3);
+ x2:=x1+ODEstep;
+ y2:=y1+ODEstep*k2;}
+ \ODEline{z1}{z2}
+ \mfsrc{
+ if ((y2>yneg) and (y2<ypos)): x1:=x2; y1:=y2 fi;
+ endfor
+ }
+ \end{connect}
+}
+%%% Integral curve using fourth order Runge--Kutta method.
+\def\trajectoryRKF#1#2{
+ \begin{connect}
+ \mfsrc{x1:=#1;y1:=#2;
+ for i=1 upto ODEstepcount:
+ k1:=ODErhs(x1,y1);
+ x3:=x1+(ODEstep/2);
+ y3:=y1+k1*(ODEstep/2);
+ k2:=ODErhs(x3,y3);
+ y4:=y1+k2*(ODEstep/2);
+ k3:=ODErhs(x3,y4);
+ y5:=y1+k3*(ODEstep/2);
+ k4:=ODErhs(x3,y5);
+ kk:=(k1+2*k2+2*k3+k4)/6;
+ x2:=x1+ODEstep;
+ y2:=y1+ODEstep*kk;}
+ \ODEline{z1}{z2}
+ \mfsrc{
+ if ((y2>yneg) and (y2<ypos)): x1:=x2; y1:=y2 fi;
+ endfor
+ }
+ \end{connect}
+}
+\def\ODEarrow#1#2{
+ \mfsrc{x1:=#1; y1:=#2;
+ x3:=x1+(ODEarrowlength)/((xscale)++(ODErhs(#1,#2)*yscale));
+ y3:=y1+(ODEarrowlength*ODErhs(#1,#2))/((xscale)++(ODErhs(#1,#2)*yscale));
+ if y3>y1:ODEcolorarrow:=blue else: ODEcolorarrow:=red fi;
+ }
+ \ifcolorODEarrow
+ \drawcolor{ODEcolorarrow} \headcolor{ODEcolorarrow}
+ \fi
+ \draw\arrow\lines{z1,z3}
+}
+
+\def\ODEarrows#1{\ODE@cycle@points#1;,;}
+\def\trajectories#1{\ODE@cycle@IC#1;,;}
+\def\ODE@last@point{}
+\def\ODE@cycle@points#1,#2;{\def\temp{#1}\ifx\temp\ODE@last@point\let\next\relax
+ \else\ODEarrow{#1}{#2}\relax\let\next\ODE@cycle@points\fi\next}
+\def\ODE@cycle@IC#1,#2;{\def\temp{#1}\ifx\temp\ODE@last@point\let\next\relax
+ \else
+ \trajectoryRKF{#1}{#2}\relax\let\next\ODE@cycle@IC\fi\next}
+\mfsrc{path p,q;color ODEcolorarrow;}
+
+%%% One-dimensional autonomous systems y'=f(y) where '=d/dx
+\def\ODEharrow#1{
+ \mfsrc{x1:=#1;
+ if ODErhs(0,x1)>0: x3:=x1+ODEarrowlength else: x3:=x1-ODEarrowlength fi;
+ if ODErhs(0,x1)*ODErhs(0,x3)<0: x1:=-100;x3:=-100 fi;
+ if x3>x1:ODEcolorarrow:=blue else: ODEcolorarrow:=red fi;
+ }
+ \ifcolorODEarrow \drawcolor{ODEcolorarrow}
+ \headcolor{ODEcolorarrow} \fi
+ \pen{1.5pt}
+ \draw\arrow\lines{(x1,0),(x3,0)}
+}
+
+\def\ODEvarrow#1{
+ \mfsrc{x1:=#1;
+ if ODErhs(0,#1)>0:
+ x3:=x1+(ODEarrowlength/yscale) else: x3:=x1-(ODEarrowlength/yscale) fi;
+ if ODErhs(0,x1)*ODErhs(0,x3)<0: x1:=-100;x3:=-100 fi;
+ if x3>x1:ODEcolorarrow:=blue else: ODEcolorarrow:=red fi;
+ }
+ \ifcolorODEarrow \drawcolor{ODEcolorarrow}
+ \headcolor{ODEcolorarrow} \fi
+ \pen{1.5pt}
+ \draw\arrow\lines{(0,x1),(0,x3)}
+}
+
+%%% Two-dimensional autonomous systems x'=f(x,y), y'=g(x,y) where '=d/dt
+\def\ASdefineequations#1#2{\fdef{ASf}{x,y}{#1}\fdef{ASg}{x,y}{#2}}
+
+\def\AStrajectory#1#2{
+ \begin{connect}
+ \mfsrc{x1:=#1;y1:=#2;
+ for i=1 upto ODEstepcount:
+ x2:=x1+ODEstep*ASf(x1,y1);
+ y2:=y1+ODEstep*ASg(x1,y1);}
+ \ODEline{z1}{z2}
+ \mfsrc{
+ if ((y2>yneg) and (y2<ypos)): x1:=x2; y1:=y2 fi;
+ endfor
+ }
+ \end{connect}
+}
+\def\ASarrow#1#2{
+ \mfsrc{x1:=#1; y1:=#2;
+ x3:=x1+(ODEarrowlength*ASf(#1,#2))/((ASf(#1,#2)*xscale)++(ASg(#1,#2)*yscale ));
+ y3:=y1+(ODEarrowlength*ASg(#1,#2))/((ASf(#1,#2)*xscale)++(ASg(#1,#2)*yscale ));
+ if y3>y1:ODEcolorarrow:=blue else: ODEcolorarrow:=red fi;
+ }
+ \ifcolorODEarrow
+ \drawcolor{ODEcolorarrow} \headcolor{ODEcolorarrow}
+ \fi
+ \draw\arrow\lines{z1,z3}
+}
+
+\def\ASarrows#1{\AS@cycle@points#1;,;}
+\def\AS@cycle@points#1,#2;{\def\temp{#1}\ifx\temp\ODE@last@point\let\next\relax
+ \else\ASarrow{#1}{#2}\relax\let\next\AS@cycle@points\fi\next}
+\def\AStrajectories#1{\AS@cycle@IC#1;,;}
+\def\AS@cycle@IC#1,#2;{\def\temp{#1}\ifx\temp\ODE@last@point\let\next\relax
+ \else
+ \AStrajectoryRKF{#1}{#2}\relax\let\next\AS@cycle@IC\fi\next}
+\def\AStrajectoryRKF#1#2{
+ \begin{connect}
+ \mfsrc{x1:=#1;y1:=#2;
+ TIMEsteps:=abs(TIMEend/TIMEstep);
+ TIME:=0;
+ for i=1 upto TIMEsteps:
+ k1:=ASf(x1,y1);
+ l1:=ASg(x1,y1);
+ k2:=ASf(x1+(TIMEstep*k1/2),y1+(TIMEstep*l1/2));
+ l2:=ASg(x1+(TIMEstep*k1/2),y1+(TIMEstep*l1/2));
+ k3:=ASf(x1+(TIMEstep*k2/2),y1+(TIMEstep*l2/2));
+ l3:=ASg(x1+(TIMEstep*k2/2),y1+(TIMEstep*l2/2));
+ k4:=ASf(x1+(TIMEstep*k3),y1+(TIMEstep*l3));
+ l4:=ASg(x1+(TIMEstep*k3),y1+(TIMEstep*l3));
+ k5:=((k1)/6)+((k2)/3)+((k3)/3)+((k4)/6);
+ l5:=(l1/6)+(l2/3)+(l3/3)+(l4/6);
+ x2:=x1+(TIMEstep*k5);
+ y2:=y1+(TIMEstep*l5);}
+ \ODEline{z1}{z2}
+ \mfsrc{
+ if ((y2>yneg) and (y2<ypos) and (x2<xpos) and (x2>xneg)): x1:=x2; y1:=y2 fi;
+ endfor
+ }
+ \end{connect}
+}
+
+\catcode`\@12\relax
+%</tex>
+%<sty>\input mfpic4ode.tex\relax
+% \end{macrocode}
+%
+% \Finale
+\endinput
+
+
+
diff --git a/macros/latex/contrib/mfpic4ode/mfpic4ode.ins b/macros/latex/contrib/mfpic4ode/mfpic4ode.ins
new file mode 100644
index 0000000000..a1fe49991f
--- /dev/null
+++ b/macros/latex/contrib/mfpic4ode/mfpic4ode.ins
@@ -0,0 +1,59 @@
+%%
+%% Copyright (C) 2007 by Robert Marik <marik@mendelu.cz>
+%%
+%% This file may be distributed and/or modified under the conditions of
+%% the LaTeX Project Public License, either version 1.2 of this license
+%% or (at your option) any later version. The latest version of this
+%% license is in:
+%%
+%% http://www.latex-project.org/lppl.txt
+%%
+%% and version 1.2 or later is part of all distributions of LaTeX version
+%% 1999/12/01 or later.
+%%
+
+\input docstrip.tex
+\keepsilent
+
+\usedir{tex/latex/mfpic4ode}
+
+\preamble
+
+This is a generated file.
+
+Copyright (C) 2007 by Robert Marik <marik@mendelu.cz>
+
+This file may be distributed and/or modified under the conditions of
+the LaTeX Project Public License, either version 1.2 of this license
+or (at your option) any later version. The latest version of this
+license is in:
+
+ http://www.latex-project.org/lppl.txt
+
+and version 1.2 or later is part of all distributions of LaTeX version
+1999/12/01 or later.
+
+\endpreamble
+
+\generate{\file{mfpic4ode.tex}{\from{mfpic4ode.dtx}{tex}}
+}
+\generate{\file{mfpic4ode.sty}{\from{mfpic4ode.dtx}{sty}}
+}
+
+\obeyspaces
+\Msg{*************************************************************}
+\Msg{* *}
+\Msg{* To finish the installation you have to move the following *}
+\Msg{* files into a directory searched by TeX: *}
+\Msg{* *}
+\Msg{* mfpic4ode.tex *}
+\Msg{* mfpic4ode.sty *}
+\Msg{* *}
+\Msg{* To produce the documentation run the file *}
+\Msg{* mfpic4ode.dtx through LaTeX. *}
+\Msg{* *}
+\Msg{* Happy TeXing! *}
+\Msg{* *}
+\Msg{*************************************************************}
+
+\endbatchfile
diff --git a/macros/latex/contrib/mfpic4ode/mfpic4ode.pdf b/macros/latex/contrib/mfpic4ode/mfpic4ode.pdf
new file mode 100644
index 0000000000..ff9a173961
--- /dev/null
+++ b/macros/latex/contrib/mfpic4ode/mfpic4ode.pdf
Binary files differ