For the first step $t_1$,
only the initial condition at $t_0$ is available.
-Therefore backward-Euler is used, i.e.,
+Therefore backward Euler is used, i.e.,
\begin{equation}\label{eq:BE}
\tfrac{1}{h_1}[q(x_1)-q(x_0)] + f(x_1) = b_1.
\end{equation}
\! + \! \frac{\alpha_2}{h_M} C_{M-2} \frac{\ud x_{M-2}}{\ud x_0}\right].
\end{equation}
+\begin{algorithm}
+\caption{The matrix-free\index{matrix-free} method for
+ Krylov subspace\index{iterative method!Krylov subspace} construction.}
+\label{alg:mf_Gear}
+ \KwIn{ current Krylov subspace basis vector $v$,
+ time step lengths $h_i$,
+ saved $C_i$ matrices and LU\index{LU} factors of $J_i$, $i=0,\ldots,M$}
+ \KwOut{ matrix-vector product $r$, such that $r = A v$}
+ solve $J_1 p_2 = h_1^{-1} C_0 v$ for $p_2$\;
+ solve $J_2 p_1 = h_2^{-1} \left( \alpha_1 C_1 p_2 + \alpha_2 C_0 v \right)$ for $p_1$\;
+ \For{$i=3 \ldots M$}{
+ solve
+ $J_i p_0 = \alpha_1 h_i^{-1} C_{i-1} p_1 + \alpha_2 h_i^{-1} C_{i-2} p_2$
+ for $p_0$ \label{line:mf_Gear_loop}\;
+ $p_2 = p_1$\;
+ $p_1 = p_0$\;
+ }
+ $r = (k-1) p_0 - k v$ \label{line:shift}\;
+\end{algorithm}
+
+
%% \begin{algorithm}
%% \caption{The matrix-free\index{matrix-free} method for
%% Krylov subspace\index{Krylov subspace} construction.}
%% \end{algorithmic}
%% \end{algorithm}
-
%It is easy to judge that, as $J$ is dense, explicitly
%computing $J$ is expensive and impractical for large problems.
%In addition, with a dense matrix,
