\% ********** Class 1D PDEs **********
\topic pde1_Class {PDEs on the interval}
The class {\bf PDEs on the interval} comprises
systems of evolutionary partial differential equations (PDEs) on the unit 
interval with separated boundary conditions that can by written in the 
following form:
\par
\qc{$u_{t}(x,t)\:=\:F(u_{xx}(x,t),\:u_{x}(x,t),\:u(x,t),\:x,\:\alpha),\:\: x \in (0,1)$}
\par
\qc{
$
\system 2 rcl 1.2 1 
{f^{0}(u_{x}(0,t),u(0,t),\alpha)}{=}{0}
{f^{1}(u_{x}(1,t),u(1,t),\alpha)}{=}{0,}
$
}
\par
where $u(x,t)\in {\f b* R}^n$ is a vector of component distributions at $(x,t)$,
$\alpha\in {\f b* R}^m$ is a vector of numerical parameters,
$F$ is a vector-function specifying the evolution of the distributions, and 
$f^{0,1}$ are boundary conditions at the left and right end point, respectively.
$F$ and $f^{0,1}$ are assumed to be sufficiently smooth but has no other special 
properties.
\par
A system of PDEs is specified or edited
in the \jump kh_RHSWin {System specification dialog box}
which is opened by the {\bf Select|System} command from the
\jump kh_MainWin {Main window.}\par
There are following class-specific fields in the dialog box.\par

{\desc 4
{\bf Components}\>is where you list the names of components.\par
{\bf Parameters}\>is where you list names of parameters your system
depends on.\par
{\bf Time}\>is where you may type in the name of time variable
which must be a simple name.\par
{\bf Space\+coordinate}\>is where you type in the name of the space variable.\par
}
\par
While specifying the system of PDEs on the interval, use the following
notations for the derivatives of component distributions:
\par
{
\desc 4
{\bf Derivative} \> {\bf Notation}\par
{\li 2 $u_{t}$} \> u\|_t\par
{\li 2 $u_{x}$} \> u\|_x\par
{\li 2 $u_{xx}$} \> u\|_x\|_x\par
}
\par
Here u is the name of a component; t is the name of time; x is the
name of the space coordinate.
\par
To specify the boundary conditions, use the corresponding command from 
the {\bf Specify} menu. While specify the boundary conditions, you have to use
the following notations in the appropriate windows:
\par
{
\desc 3
{\bf Component} \> {\bf Notation}\par
{$f^{0}_{i+1}(u(0,t),u_{x}(0,t),\alpha)=0$} \> BC\|_i=$f^{0}_{i}$(u(0),u\|_x(0),$\alpha$)\par
{$f^{1}_{i+1}(u(1,t),u_{x}(1,t),\alpha)=0$} \> BC\|_i=$f^{1}_{i}$(u(1),u\|_x(1),$\alpha$),\par
}
\par
where $i=0,1,\ldots,n-1.$
There are additional commands in the {\bf Specify} menu of the
\jump kh_RHSWin {System specification dialog box} for this class.
Use them to specify boundary conditions and their derivatives
(if you chose to supply derivatives by routines).
There is one minor difference in \jump kh_UserDer {derivatives notation:}
the first parameter to the macros ${\bf der}k$ must be the ordinal
number (starting from 0) not the name of a varivable.
\par
{\bf Example:} For the system
\par
\qc{$u_{t}(x,t)=u_{xx}(x,t)+\lambda e^{u(x,t)},\:\: x \in (0,1)$}
\par
\qc{$u(0,t)=0,\: u(1,t)=0,$}
\par
the following specification is valid.
\par
{
\desc 7
{\bf Components:} \>$U$\par
{\bf Parameters:} \>$lambda$\par
{\bf Time:}\> $t$\par
{\bf Space coordinate:} \>$x$\par
{\bf RHS:} \>U\|_t=U\|_x\|_x+lambda*exp(U);\par
{\bf Left end BC:}\> BC\|_0=U(0);\par
{\bf Right end BC:} \>BC\|_0=U(1);\par
}
\par
If you selected to specify the first-order derivatives yourself, you can write
for nonzero values
\par
{
\desc 5
{\bf RHS derivatives:} \>der1(U,U)=lambda*exp(U);\par
 \>der1(U,U\|_x\|_x)=1;\par
 \>der1(U,lambda)=exp(U);\par
{\bf Left end BC derivatives:} \>der1(0,U(0))=1;\par
{\bf Right end BC derivatives:} \>der1(0,U(1))=1;\par
}
\par
The following constants are defined and may be used
in RHS, its derivatives and local functions:\par
{\desc 2
{\bf COMPONENTSDIM}\>the number of components.\par
{\bf PARAMETERSDIM}\>the number of parameters.\par
}
\end
\% pde1_Class

\topic pde1_Types {types of curves}
CONTENT currently supports the continuation of the following curves:
\par
\par
\jump pd1_or {orbit}
\par
\jump pd1_ss {steady state curve}
\end
\% pde1_Types 


\topic pd1_integr {integrators}
The computation of orbits is developed and implemented in cooperation with
{\bf A.R. Skovoroda} (Institute of Mathematical Problems of Biology, Pushchino, Moscow Region, Russia)
and {\bf W. Govaerts} and {\bf B. Sijnave} (Vakgroep Toegepaste Wiskunde en
Informatica, Universiteit Gent, Belgium).
\par
Given a \jump pde1_Class {system of PDEs} on the unit interval, 
the {\bf orbit} is a curve parametrized by time in the space of smooth
vector functions on the unit interval, satisfying the following initial-boundary-value problem
with separated boundary conditions:
\par
\qc{
$
\system 3 rcl 1.2 1 
{u_{t}(x,t)}{=}{F(u_{xx}(x,t),\:u_{x}(x,t),\:u(x,t),\:x,\:\alpha),\:\: x \in (0,1)}
{f^{0}(u_{x}(0,t),u(0,t),\alpha)}{=}{0}
{f^{1}(u_{x}(1,t),u(1,t),\alpha)}{=}{0,}
$
}
\qr {(1)}
\par
where $F:\: {\f b* R}^{3n+1+m} \rightarrow {\f b* R}^n,\: 
f^{(0,1)}:\: {\f b* R}^{2n+m} \rightarrow {\f b* R}^n.$
\par
The initial-boundary-value problem (1) is substituted in CONTENT by a finite-difference 
approximation as following. Introduce a nonuniform mesh
\par 
\qc{
$
\Theta_N=\{0=x_0<x_1<\ldots<x_{N-1}<x_N=1\}
$
}
\qr {(2)}
\par
with the number of mesh intervals $N>6$ and set $h_i=x_{i}-x_{i-1}$ for $i=1,\ldots,N$.
Denote the corresponding mesh values of the solution $u(x_i,t)$ by $u_i(t),\ i=0,1,\ldots,N$.
\par
For the left end point we take
\par
\qc{
$
u_x(0,t)\approx {\tilde u}_x(0,t)={\hat \varphi}_{0}u_{0}(t)+{\hat \varphi}_{1}u_1(t)+{\hat \varphi}_{2}u_2(t),
$
}
\qr {(3)}
\par
where 
\par
\qc{
$
{\hat \varphi}_0=-\frac{2h_1+h_2}{h_{1}(h_1+h_2)},\: 
{\hat \varphi}_1=\frac{h_1+h_2}{h_{1}h_2},\: 
{\hat \varphi}_2=-\frac{h_1}{h_{2}(h_1+h_2)}.
$
}
\par
Associate to each internal mesh point with $i=1,2,\ldots,N-1,$
\par
\qc{
$
{\tilde x}_i=x_i+\frac{1}{3}(h_{i+1}-h_i),\ 
$
}
\par
and take 
\par
\qc{
$
u({\tilde x}_i,t)\approx{\tilde u}({\tilde x}_i,t)\ =\ \eta_{i-1}u_{i-1}(t)+\eta_{i}u_{i}(t)+\eta_{i+1}u_{i+1}(t),
$
}
\qr {(4)}
\par
\qc{
$
u_x({\tilde x}_i,t)\approx{\tilde u}_x({\tilde x}_i,t)\ =\ \varphi_{i-1}u_{i-1}(t)+\varphi_{i}u_{i}(t)
+\varphi_{i+1}u_{i+1}(t),
$
}
\qr {(5)}
\par
\qc{
$
u_{xx}({\tilde x}_i,t)\approx{\tilde u}_{xx}({\tilde x}_i,t)\ =\ \psi_{i-1}u_{i-1}(t)+\psi_{i}u_{i}(t)
+\psi_{i+1}u_{i+1}(t),
$
}
\qr {(6)}
\par
where 
\par
\qc{
$
\eta_{i-1}=-\frac{(h_{i+1}-h_i)(2h_{i+1}+h_i)}{9h_{i}(h_i+h_{i+1})},\: 
\eta_i=\frac{(2h_i+h_{i+1})(2h_{i+1}+h_i)}{9h_{i}h_{i+1}},\: 
\eta_{i+1}=\frac{(h_{i+1}-h_i)(2h_i+h_{i+1})}{9h_{i+1}(h_i+h_{i+1})},
$
}
\par
\qc{
$
\varphi_{i-1}=-\frac{2h_i+h_{i+1}}{3h_{i}(h_i+h_{i+1})},\: 
\varphi_i=\frac{h_{i+1}-h_i}{3h_{i}h_{i+1}},\: 
\varphi_{i+1}=\frac{h_i+2h_{i+1}}{3h_{i+1}(h_i+h_{i+1})},
$
}
\par
and
\par
\qc{
$
\psi_{i-1}=\frac{2}{h_{i}(h_i+h_{i+1})},\: 
\psi_i=-\frac{2}{h_{i}h_{i+1}},\: 
\psi_{i+1}=\frac{2}{h_{i+1}(h_i+h_{i+1})}.
$
}
\par
Finally, for the right end point we take
\par
\qc{
$
u_x(1,t)\approx {\tilde u}_x(1,t)={\hat \varphi}_{N-2}u_{N-2}(t)+{\hat \varphi}_{N-1}u_{N-1}(t)
+{\hat \varphi}_{N}u_N(t),
$
}
\qr {(7)}
\par
where
\par
\qc{
$
{\hat \varphi}_N=\frac{2h_N+h_{N-1}}{h_N(h_{N-1}+h_N)},\: 
{\hat \varphi}_{N-1}=-\frac{h_{N-1}+h_N}{h_{N-1}h_N},\:
{\hat \varphi}_{N-2}=\frac{h_N}{h_{N-1}(h_{N-1}+h_N)}.
$
}
\par
The above approximations have third-order accuracy for the function $u(x,t)$ and second-order 
accuracy for its first and second spatial derivatives (when $h=\max {i=1,\ldots,N}{h_i} \rightarrow \:0$).
\par
Substituting the approximations (2)-(7) into the problem (1), 
we obtain its $\Theta_N$-discretization in the form 
\par
\qc{
$
{\bf M}\dot{Y}={\bf F}(Y,\alpha)
$
}
\qr {(8)}
\par
where
\par
\qc{
$
Y=\matrix 5 l 0 0 {u_0(t)}{u_1(t)}{\ldots}{u_{N-1}(t)}{u_N(t)} ,
$
}
\par
\qc{
$
{\bf M}=\matrix 7 cccccccc 0 0
{0}{0}{0}{}{}{}{}{}
{\eta_0}{\eta_1}{\eta_2}{}{}{}{}{}
{}{\eta_1}{\eta_2}{\eta_3}{}{}{}{}
{...}{}{}{}{}{}{}{}
{}{}{}{}{\eta_{N-3}}{\eta_{N-2}}{\eta_{N-1}}{}
{}{}{}{}{}{\eta_{N-2}}{\eta_{N-1}}{\eta_{N}}
{}{}{}{}{}{0}{0}{0} ,
$
}
\par
and
\par
\qc{
$
{\bf F}(Y(t),\alpha)=\matrix 7 l 0 0 
{f^0({\tilde u}_x(0,t),u_0(t),\alpha)}
{F({\tilde u}_{xx}({\tilde x}_1,t),{\tilde u}_x({\tilde x}_1,t),{\tilde u}({\tilde x}_1,t),{\tilde x}_1,\alpha)}
{F({\tilde u}_{xx}({\tilde x}_2,t),{\tilde u}_x({\tilde x}_2,t),{\tilde u}({\tilde x}_2,t),{\tilde x}_2,\alpha)}
{\ldots}
{F({\tilde u}_{xx}({\tilde x}_{N-2},t),{\tilde u}_x({\tilde x}_{N-2},t),
{\tilde u}({\tilde x}_{N-2},t),{\tilde x}_{N-2},\alpha)}
{F({\tilde u}_{xx}({\tilde x}_{N-1},t),{\tilde u}_x({\tilde x}_{N-1},t),
{\tilde u}({\tilde x}_{N-1},t),{\tilde x}_{N-1},\alpha)}
{f^1({\tilde u}_x(1,t),u_N(t),\alpha)}
$
}
\par 
The differential-algebraic system (8) is solved by one of the following methods:
\par
\par
$\;\;\;$ \jump method_EU {Euler};
\par
$\;\;\;$ \jump method_CN {Crank-Nicolson}.
\par
\par
Since the accuracy of the above finite-difference 
approximation strongly depends on the spatial behavior of the  solution, the mesh $\Theta_N$
is {\bf adapted} along the steady state curve according to the requirement for the approximation error
for the first derivative to be uniformly distributed over the space interval. This leads to
the following mesh point selection criterium
\par
\qc{
$
S(x_{i})(x_{i}-x_{i-1})=C,
$
}
\par
where $i=1,2,\ldots,N$ and $S(x)=\sqrt{||u_{xxx}(x)||}$ is computed numerically. To recompute
the approximate solution of (1) using the new mesh, two types of interpolation are employed: 
either the interpolation by global cubic spline or the linear interpolation.
\par
Scalar functions can be {\bf monitored} along an orbit. If such a function
changes sign at a time-step, its zero can be {\bf detected} and 
computed via bi-sections to withing given tolerance.
\par
\par
An orbit can be computed form a user-supplied \jump dist_ss {distribution}. 
For this, one has to set relevant {\bf parameters to start from}
\par
\par
$\;\;\;$ \jump start_dsor {a distribution}\par
\par
\par
via the \jump kh_StarterWin {starter window}.
\par
\par
Parameters of the integrator can be modified via the corresponding
\jump kh_GeneratorWin {generator window}. \jump start_dsor {Initial data} can be
set in the corresponding \jump kh_StarterWin {starter window}.
\end	\%pde1_integr

\topic method_EU {Euler integrator}
The \jump pd1_integr {finite-difference approximation} of a \jump pde1_Class {system of PDEs}
leads to a differential-algebraic equations:
\par
\qc{
$
{\bf M}\dot{Y}(t)={\bf F}(Y(t),\alpha).
$
}
\par
Let $Y^0$ be the solution vector at time $t$. Then, an $O(h)$-approximate solution
vector $Y^1$ at time $t+h$ can be found by solving
\par
\qc{
$
{\bf M}(Y^1-Y^0)=h{\bf F}(Y^1,\alpha)
$
}
\qr {(1)}
\par
or, equivalently,
\par
\qc{
$
{\bf G}(Y^1) = {\bf F}(Y^1,\alpha) - \frac{1}{h}{\bf M}Y^1 + \frac{1}{h}{\bf M}Y^0 =0.
$
}
\qr {(2)}
\par
This is a system of nonlinear equations for $Y^1$ that we solve by Newton's method.
The Jacobian matrix of (2) has the form
\par
\qc{
$
{\bf G}_Y={\bf F}_Y - \frac{1}{h}{\bf M}.
$
}
\par
For equations (1) appearing from discretization of systems of PDEs, 
the Jacobian matrix ${\bf G}_Y$ has a special structure, therefore a {\bf modified block elimination method} 
is used to solve appearing linear equations.
\par
\jump integrpar_EU {Parameters of Euler integrator} can be modified via the corresponding
\jump kh_GeneratorWin {generator window}. \jump start_dsor {Initial data} can be
set in the corresponding \jump kh_StarterWin {starter window}.

\end

\topic method_CN {Crank-Nikolson integrator}
The \jump pd1_integr {finite-difference approximation} of a \jump pde1_Class {system of PDEs}
leads to a differential-algebraic equations:
\par
\qc{
$
{\bf M}\dot{Y}(t)={\bf F}(Y(t),\alpha)
$
}
\qr {(1)}
\par
Let $Y^0$ be the solution vector at time $t$. Then, an $O(h^2)$-approximate solution
vector $Y^1$ at time $t+h$ can be found by solving
\par
\qc{
$
{\bf M}(Y^1-Y^0)=\frac{h}{2}({\bf F}(Y^0,\alpha)+{\bf F}(Y^1,\alpha))
$
}
\par
or, equivalently,
\par
\qc{
$
{\bf G}(Y^1) = {\bf F}(Y^1,\alpha) - \frac{2}{h}{\bf M}Y^1 +
[ {\bf F}(Y^0,\alpha) + \frac{2}{h}{\bf M}Y^0 ]=0.
$
}
\qr {(2)}
\par
This is a system of nonlinear equations for $Y^1$ that we
solve by Newton's method.
The Jacobian matrix of (2) has the form
\par
\qc{
$
{\bf G}_Y={\bf F}_Y - \frac{2}{h}{\bf M}.
$
}
\par
For equation (1) appearing from discretization of systems of PDEs, 
the Jacobian matrix ${\bf G}_Y$ has a special structure, therefore a {\bf modified block elimination method} 
is used to solve appearing linear equations.
\par
\jump integrpar_CN {Parameters of Crank-Nikolson integrator} can be modified via the corresponding
\jump kh_GeneratorWin {generator window}. \jump start_dsor {Initial data} can be
set in the corresponding \jump kh_StarterWin {starter window}.

\end

\topic integrpar_EU {integrator parameters (Euler method)}
The \jump method_EU {Euler integrator} has the following parameters:
\par
\par
\desc 9
{\bf Integration data}\>\par
{\li 1 Interval} \> Integration interval length $T_{max}$\par
{\li 1 InitStepsize} \>  Initial stepsize $H_0$\par
{\li 1 Tolerance} \> Absolute step tolerance $\varepsilon$\par
{\li 1 MaxNewtonIter} \> Maximum number of Newton iterations\par
{\li 1 BetweenAdaptation} \>Number of time-steps between mesh adaptations.\par
{\bf User function zero}\>\par
{\li 1 UserFuncTolerance} \> Accuracy of zero location of user functions\par
{\li 1 MaxIterFunc} \> Maximal number of bi-sections for zero location\par
\par
The parameters can be modified via the \jump kh_GeneratorWin {generator window}.
\end

\topic integrpar_CN {integrator parameters (Crank-Nicolson method)}
The \jump method_CN {Crank-Nicolson integrator} has the following parameters:
\par
\par
\desc 11
{\bf Integration data}\>\par
{\li 1 Interval} \> Integration interval length $T_{max}$\par
{\li 1 InitStepsize} \>  Initial stepsize $H_0$\par
{\li 1 Tolerance} \> Absolute step tolerance $\varepsilon$\par
{\li 1 MaxNewtonIter} \> Maximum number of Newton iterations\par
{\li 1 MinStepsize} \> Minimal stepsize $H_{min}$\par
{\li 1 MaxStepsize} \> Maximal stepsize $H_{max}$\par
{\li 1 BetweenAdaptation} \>Number of time-steps between mesh adaptations.\par
{\bf User function zero}\>\par
{\li 1 UserFuncTolerance} \> Accuracy of zero location of user functions\par
{\li 1 MaxIterFunc} \> Maximal number of bi-sections for zero location\par
\par
The parameters can be modified via the \jump kh_GeneratorWin {generator window}.
\end

\topic pd1_ss {steady state curve}
The continuation of this curve is developed and implemented in cooperation with
A.R. Skovoroda (Institute of Mathematical Problems of Biology, Pushchino, Moscow Region, Russia).
\par
Given a \jump pde1_Class {system of PDEs} on the unit interval, 
\par
\qc{$u_{t}(x,t)\:=\:F(u_{xx}(x,t),\:u_{x}(x,t),\:u(x,t),\:x,\:\alpha),\:\: x \in (0,1),$}
\par
\qc{
$
\system 2 rcl 1.2 1 
{f^{0}(u_{x}(0,t),u(0,t),\alpha)}{=}{0,}
{f^{1}(u_{x}(1,t),u(1,t),\alpha)}{=}{0,}
$
}
\par
the {\bf steady state curve} is a one-dimensional manifold in the space of smooth
vector functions on the unit interval, satisfying the following boundary-value problem
with separated boundary conditions:
\par
\qc{
$
\system 3 rcl 1.2 1 
{F(u_{xx}(x),\:u_{x}(x),\:u(x),\:x,\:\alpha)}{=}{0,\:\: x \in (0,1)}
{f^{0}(u_{x}(0),u(0),\alpha)}{=}{0}
{f^{1}(u_{x}(1),u(1),\alpha)}{=}{0,}
$
}
\qr {(1)}
\par
where $F:\: {\f b* R}^{3n+1+m} \rightarrow {\f b* R}^n,\: 
f^{(0,1)}:\: {\f b* R}^{2n+m} \rightarrow {\f b* R}^n.$
\par
The boundary-value problem (1) is substituted in CONTENT by a finite-difference approximation as following.
Introduce a nonuniform mesh
\par 
\qc{
$
\Theta_N=\{0=x_0<x_1<\ldots<x_{N-1}<x_N=1\}
$
}
\qr {(2)}
\par
with the number of mesh intervals $N>6$ and set $h_i=x_{i}-x_{i-1}$ for $i=1,\ldots,N$.
Denote the corresponding mesh values of the solution $u(x_i)$ by $u_i,\ i=0,1,\ldots,N$.
\par
For the left end point we take
\par
\qc{
$
u_x(0)\approx {\tilde u}_x(0)={\hat \varphi}_{0}u_{0}+{\hat \varphi}_{1}u_1+{\hat \varphi}_{2}u_2,
$
}
\qr {(3)}
\par
where 
\par
\qc{
$
{\hat \varphi}_0=-\frac{2h_1+h_2}{h_{1}(h_1+h_2)},\: 
{\hat \varphi}_1=\frac{h_1+h_2}{h_{1}h_2},\: 
{\hat \varphi}_2=-\frac{h_1}{h_{2}(h_1+h_2)}.
$
}
\par
Associate to each internal mesh point with $i=1,2,\ldots,N-1,$
\par
\qc{
$
{\tilde x}_i=x_i+\frac{1}{3}(h_{i+1}-h_i),\ 
$
}
\par
and take 
\par
\qc{
$
u({\tilde x}_i)\approx{\tilde u}({\tilde x}_i)\ =\ \eta_{i-1}u_{i-1}+\eta_{i}u_{i}+\eta_{i+1}u_{i+1},
$
}
\qr {(4)}
\par
\qc{
$
u_x({\tilde x}_i)\approx{\tilde u}_x({\tilde x}_i)\ =\ \varphi_{i-1}u_{i-1}+\varphi_{i}u_{i}+\varphi_{i+1}u_{i+1},
$
}
\qr {(5)}
\par
\qc{
$
u_{xx}({\tilde x}_i)\approx{\tilde u}_{xx}({\tilde x}_i)\ =\ \psi_{i-1}u_{i-1}+\psi_{i}u_{i}+\psi_{i+1}u_{i+1},
$
}
\qr {(6)}
\par
where 
\par
\qc{
$
\eta_{i-1}=-\frac{(h_{i+1}-h_i)(2h_{i+1}+h_i)}{9h_{i}(h_i+h_{i+1})},\: 
\eta_i=\frac{(2h_i+h_{i+1})(2h_{i+1}+h_i)}{9h_{i}h_{i+1}},\: 
\eta_{i+1}=\frac{(h_{i+1}-h_i)(2h_i+h_{i+1})}{9h_{i+1}(h_i+h_{i+1})},
$
}
\par
\qc{
$
\varphi_{i-1}=-\frac{2h_i+h_{i+1}}{3h_{i}(h_i+h_{i+1})},\: 
\varphi_i=\frac{h_{i+1}-h_i}{3h_{i}h_{i+1}},\: 
\varphi_{i+1}=\frac{h_i+2h_{i+1}}{3h_{i+1}(h_i+h_{i+1})},
$
}
\par
and
\par
\qc{
$
\psi_{i-1}=\frac{2}{h_{i}(h_i+h_{i+1})},\: 
\psi_i=-\frac{2}{h_{i}h_{i+1}},\: 
\psi_{i+1}=\frac{2}{h_{i+1}(h_i+h_{i+1})}.
$
}
\par
Finally, for the right end point we take
\par
\qc{
$
u_x(1)\approx {\tilde u}_x(1)={\hat \varphi}_{N-2}u_{N-2}+{\hat \varphi}_{N-1}u_{N-1}+{\hat \varphi}_{N}u_N,
$
}
\qr {(7)}
\par
where
\par
\qc{
$
{\hat \varphi}_N=\frac{2h_N+h_{N-1}}{h_N(h_{N-1}+h_N)},\: 
{\hat \varphi}_{N-1}=-\frac{h_{N-1}+h_N}{h_{N-1}h_N},\:
{\hat \varphi}_{N-2}=\frac{h_N}{h_{N-1}(h_{N-1}+h_N)}.
$
}
\par
The above approximations have third-order accuracy for the function $u(x)$ and second-order 
accuracy for its first and second derivatives (when $h=\max {i=1,\ldots,N}{h_i} \rightarrow \:0$).
\par
Substituting the approximations (2)-(7) into the problem (1), 
we obtain its $\Theta_N$-discretization in the form 
\par
\qc{
$
{\bf F}(Y)=0
$
}
\qr {(8)}
\par
where
\par
\qc{
$
Y=\matrix 6 l 0 0 {u_0}{u_1}{\ldots}{u_{N-1}}{u_N}{\alpha}
$
}
\par
and
\par
\qc{
$
{\bf F}(Y)=\matrix 7 l 0 0 
{f^0({\tilde u}_x(0),u_0,\alpha)}
{F({\tilde u}_{xx}({\tilde x}_1),{\tilde u}_x({\tilde x}_1),{\tilde u}({\tilde x}_1),{\tilde x}_1,\alpha)}
{F({\tilde u}_{xx}({\tilde x}_2),{\tilde u}_x({\tilde x}_2),{\tilde u}({\tilde x}_2),{\tilde x}_2,\alpha)}
{\ldots}
{F({\tilde u}_{xx}({\tilde x}_{N-2}),{\tilde u}_x({\tilde x}_{N-2}),
{\tilde u}({\tilde x}_{N-2}),{\tilde x}_{N-2},\alpha)}
{F({\tilde u}_{xx}({\tilde x}_{N-1}),{\tilde u}_x({\tilde x}_{N-1}),
{\tilde u}({\tilde x}_{N-1}),{\tilde x}_{N-1},\alpha)}
{f^1({\tilde u}_x(1),u_N,\alpha)}
$
}
\par 
The system (8) specifies a curve in ${\bf R}^K$, with $K = \dim Y = (N+1)n+m$. The components of ${\bf F}$
in (8) are the \jump kh_DefFunc {defining functions} of the {\bf steady state curve}.
The Jacobian matrix ${\bf F}_Y$ of (8) involved in the continuation
has a special structure, therefore a {\bf modified block elimination method} 
is used to solve appearing linear equations.
\par
Since the accuracy of the above finite-difference 
approximation strongly depends on the spatial behavior of the steady-state solution, the mesh $\Theta_N$
is {\bf adapted} along the steady state curve according to the requirement for the approximation error
for the first derivative to be uniformly distributed over the space interval. This leads to
the following mesh point selection criterium
\par
\qc{
$
S(x_{i})(x_{i}-x_{i-1})=C,
$
}
\par
where $i=1,2,\ldots,N$ and $S(x)=\sqrt{||u_{xxx}(x)||}$ is computed numerically. To recompute
the approximate solution of (1) using the new mesh, two types of interpolation are employed: 
either the interpolation by global cubic spline or the linear interpolation.
\par
Several \jump test_ss {test functions} can be computed along the curve to
detect and process \jump test_ss {steady state singularities}.
\par
\par
A steady state curve can be continued form a user-supplied \jump dist_ss {distribution} or 
a \jump ss_ss {steady state}. To \jump conti {continue} the steady state curve, one needs to set relevant
{\bf parameters to start from}
\par
\par
$\;\;\;$ \jump start_dss {a distribution}\par
$\;\;\;$ \jump start_eqe {a steady state}\par
\par
\par
via the \jump kh_StarterWin {starter window}.
\end	\%pde1_ss

\topic dist_ss {distribution}
A vector function $u(x)$ on the unit interval $[0,1]$ is called the {\bf distribution}.
In CONTENT the distribution $u(x)$ is approximated by its values at mesh points
\par
\qc{
$
\Theta_N=\{0=x_0<x_1<\ldots<x_{N-1}<x_N=1\}.
$
}
\par
You can compute 
\par
{\li 2
\desc 2
\jump pd1_or {an orbit}\>
 \par
\jump pd1_ss {a steady state curve}\>
 \par
}
\par
starting from a distribution at some parameter values. In the last case, 
the distribution should be rather close to a steady state.
\end

\topic pd1_or {orbit}
Given a \jump pde1_Class {system of PDEs} on the unit interval, 
\par
\qc{$u_{t}(x,t)\:=\:F(u_{xx}(x,t),\:u_{x}(x,t),\:u(x,t),\:x,\:\alpha),\:\: x \in (0,1),$}
\par
\qc{
$
\system 2 rcl 1.2 1 
{f^{0}(u_{x}(0,t),u(0,t),\alpha)}{=}{0,}
{f^{1}(u_{x}(1,t),u(1,t),\alpha)}{=}{0,}
$
}
\par
an {\bf orbit} is defined by the solution $u=u(x,t)$ as a parametrized by 
time curve starting at $u(x,t_0)=u_0(x)$ and belonging to the function space 
of smooth vector functions on the unit interval $[0,1]$.
\par
In CONTENT orbits are generated numerically by \jump pd1_integr {integrator}.
\end

\topic ss_ss {steady state}
Given a \jump pde1_Class {system of PDEs} on the unit interval, 
\par
\qc{$u_{t}(x,t)\:=\:F(u_{xx}(x,t),\:u_{x}(x,t),\:u(x,t),\:x,\:\alpha),\:\: x \in (0,1),$}
\par
\qc{
$
\system 2 rcl 1.2 1 
{f^{0}(u_{x}(0,t),u(0,t),\alpha)}{=}{0,}
{f^{1}(u_{x}(1,t),u(1,t),\alpha)}{=}{0,}
$
}
\par
a {\bf steady state} is characterized by a
vector function $u(x)$ on the unit interval, satisfying the following boundary-value problem:
\par
\qc{
$
\system 3 rcl 1.2 1 
{F(u_{xx}(x),\:u_{x}(x),\:u(x),\:x,\:\alpha)}{=}{0,\:\: x \in (0,1)}
{f^{0}(u_{x}(0),u(0),\alpha)}{=}{0}
{f^{1}(u_{x}(1),u(1),\alpha)}{=}{0,}
$
}
\par
where $F:\: {\f b* R}^{3n+1+m} \rightarrow {\f b* R}^n,\: 
f^{(0,1)}:\: {\f b* R}^{2n+m} \rightarrow {\f b* R}^n.$
\par
In CONTENT the steady state $u(x)$ is approximated by its values at mesh points
\par
\qc{
$
\Theta_N=\{0=x_0<x_1<\ldots<x_{N-1}<x_N=1\}.
$
}
\end

\topic test_ss {test functions and singularities along steady state curve}
The following \jump kh_TestFunc {test-functions} can be computed along the \jump pd1_ss {steady state curve}
${\bf F}(Y)=0$:
\par
\qc{
$
\table 2 rcl 0 0 
{\psi_1}{\:=\:}{det \matrix 2 c 0 0 {{\bf F}_Y}{v^T},}
{\psi_2}{\:=\:}{v_{n+1}.}
$}
\par
Here $v=(v_1,v_2,\ldots,v_{K})\in {\f b* R}^{K}$ is the normalized tangent vector to the curve
at a point $Y$ on it.
\par
\par
The following {\f b* singularities} can be detected and located as regular zeroes of
the steady state test-functions:
\par
{\li 2
\desc 2 
{\bf Branching point}:\> $\psi_1=0.$
\par
{\bf Limit point}:\> $\psi_2=0,\; \psi_1 \neq 0.$
\par
}
\end \% 

\topic start_dss {starter parameters (steady state from distribution)}
The \jump kh_Starter {starter} which produces necessary data to \jump kh_Generator {generate} 
the \jump pd1_ss {steady state curve} from a user-supplied {\bf distribution} has the following parameters
that can be modified via the \jump kh_StarterWin {starter window.}
\par
\par
{
\desc 22
{\bf Initial point}\>\par
{\li 1 $t$} \> A name of time variable and its value at the initial point.\par
{\li 1 $\alpha$} \> A list of parameter names and their values at the initial point. 
 The names of active parameters are highlighted.\par
{\bf Discretization data}\> \par
{\li 1 MeshPoints} \> The number of mesh points.\par
{\bf Corrector data}\>\par
{\li 1 MaxIter} \>  Maximal number of corrections to locate the first point.\par
{\li 1 MaxNewtonIter} \>  Maximal number of corrections with recomputation of the Jacobian
                        matrix.\par
{\li 1 VarTolerance} \>  Tolerance with respect to phase variables and parameters.\par
{\li 1 FuncTolerance} \>  Tolerance with respect to defining functions.\par
{\bf Jacobian data}\> \par
{\li 1 Increment} \> Increment to approximate partial derivatives by finite differences.\par
{\bf Monitor singularities}\>\par
{\li 1 \jump test_ss {Singularity types}}\> Monitor toggle ({\bf yes}|{\bf no}).\par
{\bf Standard functions}\>\par
{\li 1 \jump pd1_stdf {Function name}}\> Compute toggle ({\bf yes}|{\bf no}).\par
{\bf User defined functions}\>\par
{\li 1 \jump kh_UserFunc {Function name}}\> Process toggle ({\bf ignore}|{\bf monitor}|{\bf detect}|{\bf append}).\par
{\bf Interpolation method}\> \par
{\li 1 Interpolation} \>   Method toggle ({\bf linear}|{\bf cubic_spline}).\par 
{\bf Set initial point}\> \par
{\li 1 SetInitPoint} \> \jump kh_FuncPar {Functional parameter} assigning initial steady state
 and parameters. When it is specified and
 activated, it is called by the starter before computing the first point.  Otherwise, CONTENT
 uses the values corresponding to the current point.  \par
}
\par
\include pd1_init
\end

\topic start_dsor {starter parameters (orbit from distribution)}
The \jump kh_Starter {starter} which produces necessary data to \jump kh_Generator {generate} 
the \jump pd1_or {orbit} from a \jump dist_ss {distribution} has the following parameters
that can be modified via the \jump kh_StarterWin {starter window.}
\par
\par
{
\desc 15
{\bf Initial point}\>\par
{\li 1 $t$} \> A name of time variable and its value at the initial point.\par
{\li 1 $\alpha$} \> A list of parameter names and their values at the initial point.\par
{\bf Discretization data}\> \par
{\li 1 MeshPoints} \> The number of mesh points.\par
{\bf Jacobian data}\> \par
{\li 1 Increment} \> Increment to approximate partial derivatives by finite differences.\par
{\bf Standard functions}\>\par
{\li 1 \jump pd1_stdf {Function name}}\> Compute toggle ({\bf yes}|{\bf no}).\par
{\bf User defined functions}\>\par
{\li 1 \jump kh_UserFunc {Function name}}\> Process toggle ({\bf ignore}|{\bf monitor}|{\bf detect}|{\bf append}).\par
{\bf Interpolation method}\> \par
{\li 1 Interpolation} \>   Method toggle ({\bf linear}|{\bf cubic_spline}).\par 
{\bf Set initial point}\> \par
{\li 1 SetInitPoint} \> \jump kh_FuncPar {Functional parameter} assigning initial 
 distribution and parameters. When it is specified and
 activated, it is called by the starter before computing the first point.  Otherwise, CONTENT
 uses the values corresponding to the current point.  \par
}
\par
\include pd1_init
\end

\topic start_ssss {starter parameters (steady state from steady state)}
The \jump kh_Starter {starter} which produces necessary data to \jump kh_Generator {generate} 
the \jump pd1_ss {steady state curve} from a {\bf steady state} has the following parameters
that can be modified via the \jump kh_StarterWin {starter window.}
\par
\par
{
\desc 17
{\bf Initial point}\>\par
{\li 1 $t$} \> A name of time variable and its value at the initial point.\par
{\li 1 $\alpha$} \> A list of parameter names and their values at the initial point. 
 The names of active parameters are highlighted.\par
{\bf Discretization data}\> \par
{\li 1 MeshPoints} \> The number of mesh points.\par
{\bf Jacobian data}\> \par
{\li 1 Increment} \> Increment to approximate partial derivatives by finite differences.\par
{\bf Monitor singularities}\>\par
{\li 1 \jump test_ss {Singularity types}}\> Monitor toggle ({\bf yes}|{\bf no}).\par
{\bf Standard functions}\>\par
{\li 1 \jump pd1_stdf {Function name}}\> Compute toggle ({\bf yes}|{\bf no}).\par
{\bf User defined functions}\>\par
{\li 1 \jump kh_UserFunc {Function name}}\> Process toggle ({\bf ignore}|{\bf monitor}|{\bf detect}|{\bf append}).\par
{\bf Interpolation method}\> \par
{\li 1 Interpolation} \>   Method toggle ({\bf linear}|{\bf cubic_spline}).\par 
{\bf Set initial point}\> \par
{\li 1 SetInitPoint} \> \jump kh_FuncPar {Functional parameter} assigning initial 
 steady state and parameters. When it is specified and
 activated, it is called by the starter before computing the first point.  Otherwise, CONTENT
 uses the values corresponding to the current point.  \par
}
\par
\include pd1_init
\end

\topic pd1_init {-}
The function SetInitPoint is called by CONTENT several times. It has one input parameter
{\bf ip} which specifies the purpose of the call: 
\par
{\desc 4
{\bf ip}\>{\bf Meaning}\par
-1\>initialize;\par
n\>assign space and component values corresponding to nth mesh point (n=0,1,...).\par  
-2\>terminate;\par}
\par
This function should return 0 if it be called more times. Return value 1 indicates 
the last mesh point.  The names of the space variable, components, and parameters 
should be used as specified in the \jump kh_RHSWin {system specification dialog box.}
\par
In the following examples, X denotes the space variable, U and V are the names of
the components, and P is the parameter.
\par 
The first example shows how to setup an initial point given by the explicit formulas:\par
\qc{
$
U(X)=P*exp(X),\:\:V(X)=\frac{1}{X+P}.
$
}
\par
on the uniform mesh. The corresponding function SetInitPoint can be specified as following:
\par
\verbatim 8
    if (ip<0) return 0;
    if (ip>=0 && ip<MeshPoints) {
      X=ip*1.0/(MeshPoints-1);
      U=P*exp(X);
      V(X)=1.0/(X+P);
      return 0;
    }
    return 1;
\par
Note that MeshPoints has the value set in the starter window. 
\par
In the second example SetInitPoint reads the initial data from a file:
\par
\verbatim 14
    static FILE *stream;
    switch (ip) {
      case -1:
        stream=fopen("init.dat","rt");
        return 0;
      case -2:
        fclose(stream);
        return 0;
      default:
        if (fscanf(stream,"%lg %lg %lg",&X,&U,&V)==3) 
          return 0;
        else 
          return 1;
    }
\par
It is assumed that the file init.dat has three columns with X,U,V values.
If the number of the lines read differs from the specified MeshPoints, CONTENT
takes the uniform mesh with the given number of mesh points and computes new values of
the initial distribution by interpolation.
\end

\topic pd1_stdf {PDE standard functions}
Given a vector-function $u(x)$ on the interval $[0,1]$, the following {\bf standard functions}
can be computed in CONTENT: $L_2$\+-norm:
\par
\qc{
$
||u||_{2}=\sqrt{\:\Int {0}{1} {(u^{2}_1(x)+u^{2}_2(x)+\ldots+u^{2}_n(x))\: dx}};
$
}
\par
minimum and maximum of the components:
\par
\qc{
$
{Min{_u}}_i=\min {x\in[0,1]}{u_i(x)},\: \: i=1,2,\ldots,n;
$
}
\par
\qc{
$
{Max{_u}}_i=\max {x\in[0,1]}{u_i(x)},\: \: i=1,2,\ldots,n.
$
}
\par
Since in CONTENT the function $u(x)$ is approximated by its values at the mesh points
\par
\qc{
$
\Theta_N=\{0=x_0<x_1<\ldots<x_{N-1}<x_N=1\},
$
}
\par
linear or cubic spline interpolation is used to compute the above functions.
\end