MA342 Handout:             Solve system of ODEs Using Octave function LSODE
  ===========================================================================

  [X, ISTATE, MSG] = lsode (FCN, X_0, T, T_CRIT)

     Solve the set of differential equations
          dx
          -- = f(x, t)
          dt
     with
          x(t_0) = x_0

     The solution is returned in the matrix X, with each row corresponding
     to an element of the vector T.  The first element of T should be t_0
     and should correspond to the initial state of the system X_0, so that
     the first row of the output is X_0.

     The first argument, FCN, is a string, inline, or function handle that
     names the function f to call to compute the vector of right hand sides
     for the set of equations.  The function must have the form

          XDOT = f (X, T)

     in which XDOT and X are vectors and T is a scalar.

     If FCN is a two-element string array or a two-element cell array of
     strings, inline functions, or function handles, the first element
     names the function f described above, and the second element names
     a function to compute the Jacobian of f.  The Jacobian function
     must have the form

          JAC = j (X, T)

     in which JAC is the matrix of partial derivatives

                      | df_1  df_1       df_1 |
                      | ----  ----  ...  ---- |
                      | dx_1  dx_2       dx_N |
                      |                       |
                      | df_2  df_2       df_2 |
                      | ----  ----  ...  ---- |
               df_i   | dx_1  dx_2       dx_N |
         jac = ---- = |                       |
               dx_j   |  .    .     .    .    |
                      |  .    .      .   .    |
                      |  .    .       .  .    |
                      |                       |
                      | df_N  df_N       df_N |
                      | ----  ----  ...  ---- |
                      | dx_1  dx_2       dx_N |

     The fourth argument is optional, and may be used to specify a set
     of times that the ODE solver should not integrate past.  It is
     useful for avoiding difficulties with singularities and points
     where there is a discontinuity in the derivative.

     After a successful computation, the value of ISTATE will be 2
     (consistent with the Fortran version of LSODE).

     If the computation is not successful, ISTATE will be something
     other than 2 and MSG will contain additional information.

     You can use the function `lsode_options' to set optional
     parameters for `lsode'.


  Examples
  ===========================================================================

  
  1. Solve the following system on [0,5] using Octave.

         x' = y * sin(t)
         y' = x + cos(t)
         x(0) = 1
         y(0) = 0

     >> f = @(x,t) [x(2)*sin(t), x(1)*cos(t)]';
     >> x0 = [1, 0]';
     >> t = linspace(0, 5, 501);
     >> x = lsode(f, x0, t);
     >> plot(t, x(:,1), 'r-', t, x(:,2), 'b-' )
     >> legend('x', 'y', 'location', 'southwest')
     >> xlabel('t')
     >> title('Solution of System')

  2. Solve

        y''' + 2y'' - y' - 2y = exp(t)    0<=t<=3

     subject to

        y(0) = 1, y'(0) = 2, y''(0) = 0.

     >> f = @(x,t) [x(2), x(3), exp(t) + dot([2,1,-2],x)];
     >> t = linspace(0, 3, 31)';
     >> x0 = [1,2,0]';
     >> x = lsode(f, x0, t);
     >> plot(t, x(:,1))
     >> xlabel('t')
     >> ylabel('y')
     >> title('Solution')