unwind_protect

  opso = page_screen_output (0);

  tp = odetest ();

  p1 = tp.p1;
  x01 = p1.x0 (100);
  t1 = linspace (0, 20, 51);

  p2 = tp.p2;
  x02 = p2.x0 (1000);
  t2 = linspace (0, 5, 51);

  printf ("Plotted results of each test problem should always be equal, but elapsed times should differ.\n\n");

  printf ("First problem: system of 100 ODEs, Jacobian has only diagonal and one subdiagonal.\n\n");

  printf ("Full Jacobian approximated with finite differences.\n");

  lsode_options ("jacobian type", "full");

  tic; [cx, state] = lsode (p1.ode, x01, t1); toc

  if (state == 2)
    printf ("lsode reports success (state == 2)\n");
  else
    error ("lsode failed")
  endif

  printf ("plotting result to figure 1\n\n");

  figure (1);
  plot (cx);

  printf ("Fully used user-supplied analytical Jacobian.\n");
  tic; [cx, state] = lsode ({p1.ode, p1.jac}, x01, t1); toc

  if (state == 2)
    printf ("lsode reports success (state == 2)\n");
  else
    error ("lsode failed")
  endif

  printf ("plotting result to figure 2\n\n");

  figure (2);
  plot (cx);

  printf ("Banded Jacobian approximated with finite differences.\n");

  lsode_options ("jacobian type", "banded");
  lsode_options ("lower jacobian subdiagonals", 1);
  lsode_options ("upper jacobian subdiagonals", 0);

  tic;[cx, state] = lsode (p1.ode, x01, t1); toc

  if (state == 2)
    printf ("lsode reports success (state == 2)\n");
  else
    error ("lsode failed")
  endif

  printf ("plotting result to figure 3\n\n");

  figure (3);
  plot (cx);

  printf ("User-supplied analytical Jacobian used in a banded way.\n");
  tic;[cx, state] = lsode ({p1.ode, p1.jac}, x01, t1); toc

  if (state == 2)
    printf ("lsode reports success (state == 2)\n");
  else
    error ("lsode failed")
  endif

  printf ("plotting result to figure 4\n\n");

  figure (4);
  plot (cx);

  printf ("Second problem: system of 1000 independent (identical) ODEs.\n\n");

  printf ("Full Jacobian approximated with finite differences.\n");

  lsode_options ("jacobian type", "full");

  tic; [cx, state] = lsode (p2.ode, x02, t2); toc

  if (state == 2)
    printf ("lsode reports success (state == 2)\n");
  else
    error ("lsode failed")
  endif

  printf ("plotting result to figure 5\n\n");

  figure (5);
  plot (cx);

  printf ("Diagonal Jacobian approximated with finite differences.\n");

  lsode_options ("jacobian type", "diagonal");

  tic; [cx, state] = lsode (p2.ode, x02, t2); toc

  if (state == 2)
    printf ("lsode reports success (state == 2)\n");
  else
    error ("lsode failed")
  endif

  printf ("plotting result to figure 6\n\n");

  figure (6);
  plot (cx);


unwind_protect_cleanup

  page_screen_output (opso);

  figure (1);

end_unwind_protect
