--- a/liboctave/DASSL.cc
+++ b/liboctave/DASSL.cc
@@ -63,12 +63,10 @@
   liw = 0;
   lrw = 0;
 
-  info  = new int [15];
-  iwork = (int *) 0;
-  rwork = (double *) 0;
+  info.resize (15);
 
   for (int i = 0; i < 15; i++)
-    info [i] = 0;
+    info.elem (i) = 0;
 }
 
 DASSL::DASSL (const ColumnVector& state, double time, DAEFunc& f)
@@ -82,12 +80,10 @@
   liw = 20 + n;
   lrw = 40 + 9*n + n*n;
 
-  info  = new int [15];
-  iwork = new int [liw];
-  rwork = new double [lrw];
+  info.resize (15);
 
   for (int i = 0; i < 15; i++)
-    info [i] = 0;
+    info.elem (i) = 0;
 }
 
 DASSL::DASSL (const ColumnVector& state, const ColumnVector& deriv,
@@ -105,19 +101,10 @@
   liw = 20 + n;
   lrw = 40 + 9*n + n*n;
 
-  info  = new int [15];
-  iwork = new int [liw];
-  rwork = new double [lrw];
+  info.resize (15);
 
   for (int i = 0; i < 15; i++)
-    info [i] = 0;
-}
-
-DASSL::~DASSL (void)
-{
-  delete [] info;
-  delete [] rwork;
-  delete [] iwork;
+    info.elem (i) = 0;
 }
 
 void
@@ -199,12 +186,26 @@
 ColumnVector
 DASSL::do_integrate (double tout)
 {
+  ColumnVector retval;
+
+  if (restart)
+    {
+      restart = 0;
+      info.elem (0) = 0;
+    }
+
+  if (iwork.length () != liw)
+    iwork.resize (liw);
+
+  if (rwork.length () != lrw)
+    rwork.resize (lrw);
+
   integration_error = 0;
 
   if (DAEFunc::jacobian_function ())
-    iwork [4] = 1;
+    iwork.elem (4) = 1;
   else
-    iwork [4] = 0;
+    iwork.elem (4) = 0;
 
   double *px    = x.fortran_vec ();
   double *pxdot = xdot.fortran_vec ();
@@ -215,93 +216,91 @@
 
   if (stop_time_set)
     {
-      info [3] = 1;
-      rwork [0] = stop_time;
+      info.elem (3) = 1;
+      rwork.elem (0) = stop_time;
     }
   else
-    info [3] = 0;
+    info.elem (3) = 0;
 
   double abs_tol = absolute_tolerance ();
   double rel_tol = relative_tolerance ();
 
   if (initial_step_size () >= 0.0)
     {
-      rwork[2] = initial_step_size ();
-      info[7] = 1;
+      rwork.elem (2) = initial_step_size ();
+      info.elem (7) = 1;
     }
   else
-    info[7] = 0;
+    info.elem (7) = 0;
 
   if (maximum_step_size () >= 0.0)
     {
-      rwork[2] = maximum_step_size ();
-      info[6] = 1;
+      rwork.elem (2) = maximum_step_size ();
+      info.elem (6) = 1;
     }
   else
-    info[6] = 0;
+    info.elem (6) = 0;
 
   double *dummy = 0;
   int *idummy = 0;
 
-  if (restart)
-    {
-      restart = 0;
-      info[0] = 0;
-    }
+  int *pinfo = info.fortran_vec ();
+  int *piwork = iwork.fortran_vec ();
+  double *prwork = rwork.fortran_vec ();
 
 // again:
 
-  F77_FCN (ddassl, DDASSL) (ddassl_f, n, t, px, pxdot, tout, info,
-			    rel_tol, abs_tol, idid, rwork, lrw, iwork,
-			    liw, dummy, idummy, ddassl_j);
+  F77_XFCN (ddassl, DDASSL, (ddassl_f, n, t, px, pxdot, tout, pinfo,
+			     rel_tol, abs_tol, idid, prwork, lrw,
+			     piwork, liw, dummy, idummy, ddassl_j));
 
-  switch (idid)
+  if (f77_exception_encountered)
+    (*current_liboctave_error_handler) ("unrecoverable error in dassl");
+  else
     {
-    case 1: // A step was successfully taken in the
-	    // intermediate-output mode. The code has not yet reached
-	    // TOUT.
-      break;
+      switch (idid)
+	{
+	case 1: // A step was successfully taken in intermediate-output
+	        // mode. The code has not yet reached TOUT.
+	case 2: // The integration to TSTOP was successfully completed
+	        // (T=TSTOP) by stepping exactly to TSTOP.
+	case 3: // The integration to TOUT was successfully completed
+	        // (T=TOUT) by stepping past TOUT.  Y(*) is obtained by
+	        // interpolation.  YPRIME(*) is obtained by interpolation.
 
-    case 2: // The integration to TSTOP was successfully completed
-	    // (T=TSTOP) by stepping exactly to TSTOP.
-      break;
-
-    case 3: // The integration to TOUT was successfully completed
-	    // (T=TOUT) by stepping past TOUT.  Y(*) is obtained by
-	    // interpolation.  YPRIME(*) is obtained by interpolation.
-      break;
+	  retval = x;
+	  t = tout;
+	  break;
 
-    case -1: // A large amount of work has been expended.  (About 500 steps).
-    case -2: // The error tolerances are too stringent.
-    case -3: // The local error test cannot be satisfied because you
-	     // specified a zero component in ATOL and the
-	     // corresponding computed solution component is zero.
-	     // Thus, a pure relative error test is impossible for
-	     // this component.
-    case -6: // DDASSL had repeated error test failures on the last
-	     // attempted step.
-    case -7: // The corrector could not converge.
-    case -8: // The matrix of partial derivatives is singular.
-    case -9: // The corrector could not converge.  There were repeated
-	     // error test failures in this step.
-    case -10: // The corrector could not converge because IRES was
-	      // equal to minus one.
-    case -11: // IRES equal to -2 was encountered and control is being
-	      // returned to the calling program.
-    case -12: // DDASSL failed to compute the initial YPRIME.
-    case -33: // The code has encountered trouble from which it cannot
-	      // recover. A message is printed explaining the trouble
-	      // and control is returned to the calling program. For
-	      // example, this occurs when invalid input is detected.
-
-    default:
-      integration_error = 1;
-      break;
+	case -1: // A large amount of work has been expended.  (~500 steps).
+	case -2: // The error tolerances are too stringent.
+	case -3: // The local error test cannot be satisfied because you
+	         // specified a zero component in ATOL and the
+		 // corresponding computed solution component is zero.
+		 // Thus, a pure relative error test is impossible for
+		 // this component.
+	case -6: // DDASSL had repeated error test failures on the last
+		 // attempted step.
+	case -7: // The corrector could not converge.
+	case -8: // The matrix of partial derivatives is singular.
+	case -9: // The corrector could not converge.  There were repeated
+		 // error test failures in this step.
+	case -10: // The corrector could not converge because IRES was
+		  // equal to minus one.
+	case -11: // IRES equal to -2 was encountered and control is being
+		  // returned to the calling program.
+	case -12: // DDASSL failed to compute the initial YPRIME.
+	case -33: // The code has encountered trouble from which it cannot
+		  // recover. A message is printed explaining the trouble
+		  // and control is returned to the calling program. For
+		  // example, this occurs when invalid input is detected.
+	default:
+	  integration_error = 1;
+	  break;
+	}
     }
 
-  t = tout;
-
-  return x;
+  return retval;
 }
 
 Matrix

--- a/liboctave/DASSL.h
+++ b/liboctave/DASSL.h
@@ -119,7 +119,7 @@
   DASSL (const ColumnVector& x, const ColumnVector& xdot,
 	 double time, DAEFunc& f);
 
- ~DASSL (void);
+  ~DASSL (void) { }
 
   void force_restart (void);
 
@@ -146,9 +146,9 @@
   int liw;  
   int lrw;
   int idid;
-  int *info;
-  int *iwork;
-  double *rwork;
+  Array<int> info;
+  Array<int> iwork;
+  Array<double> rwork;
 
   friend int ddassl_j (double *time, double *state, double *deriv,
 		       double *pd, double *cj, double *rpar, int *ipar);

--- a/liboctave/LSODE.cc
+++ b/liboctave/LSODE.cc
@@ -72,14 +72,6 @@
 
   liw = 20 + n;
   lrw = 22 + n * (9 + n);
-
-  iwork = new int [liw];
-  rwork = new double [lrw];
-  for (int i = 4; i < 9; i++)
-    {
-      iwork[i] = 0;
-      rwork[i] = 0.0;
-    }
 }
 
 LSODE::LSODE (const ColumnVector& state, double time, const ODEFunc& f)
@@ -100,20 +92,6 @@
 
   liw = 20 + n;
   lrw = 22 + n * (9 + n);
-
-  iwork = new int [liw];
-  rwork = new double [lrw];
-  for (int i = 4; i < 9; i++)
-    {
-      iwork[i] = 0;
-      rwork[i] = 0.0;
-    }
-}
-
-LSODE::~LSODE (void)
-{
-  delete [] rwork;
-  delete [] iwork;
 }
 
 void
@@ -180,6 +158,30 @@
 ColumnVector
 LSODE::do_integrate (double tout)
 {
+  ColumnVector retval;
+
+  if (restart)
+    {
+      restart = 0;
+      istate = 1;
+    }
+
+  if (iwork.length () != liw)
+    {
+      iwork.resize (liw);
+
+      for (int i = 4; i < 9; i++)
+	iwork.elem (i) = 0;
+    }
+
+  if (rwork.length () != lrw)
+    {
+      rwork.resize (lrw);
+
+      for (int i = 4; i < 9; i++)
+	rwork.elem (i) = 0;
+    }
+
   if (jac)
     method_flag = 21;
   else
@@ -199,13 +201,12 @@
 
   // Try 5000 steps before giving up.
 
-  iwork[5] = 5000;
-  int working_too_hard = 0;
+  iwork.elem (5) = 5000;
 
   if (stop_time_set)
     {
       itask = 4;
-      rwork [0] = stop_time;
+      rwork.elem (0) = stop_time;
     }
   else
     {
@@ -215,64 +216,67 @@
   double abs_tol = absolute_tolerance ();
   double rel_tol = relative_tolerance ();
 
-  rwork[4] = (initial_step_size () >= 0.0) ? initial_step_size () : 0.0;
-  rwork[5] = (maximum_step_size () >= 0.0) ? maximum_step_size () : 0.0;
-  rwork[6] = (minimum_step_size () >= 0.0) ? minimum_step_size () : 0.0;
+  rwork.elem (4) = (initial_step_size () >= 0.0) ? initial_step_size () : 0.0;
+  rwork.elem (5) = (maximum_step_size () >= 0.0) ? maximum_step_size () : 0.0;
+  rwork.elem (6) = (minimum_step_size () >= 0.0) ? minimum_step_size () : 0.0;
 
-  if (restart)
-    {
-      restart = 0;
-      istate = 1;
-    }
+  int *piwork = iwork.fortran_vec ();
+  double *prwork = rwork.fortran_vec ();
+
+  working_too_hard = 0;
 
  again:
 
-  F77_FCN (lsode, LSODE) (lsode_f, n, xp, t, tout, itol, rel_tol,
-			  abs_tol, itask, istate, iopt, rwork, lrw,
-			  iwork, liw, lsode_j, method_flag);
+  F77_XFCN (lsode, LSODE, (lsode_f, n, xp, t, tout, itol, rel_tol,
+			   abs_tol, itask, istate, iopt, prwork, lrw,
+			   piwork, liw, lsode_j, method_flag));
 
-  switch (istate)
+  if (f77_exception_encountered)
+    (*current_liboctave_error_handler) ("unrecoverable error in lsode");
+  else
     {
-    case -13: // Return requested in user-supplied function.
-    case -6:  // error weight became zero during problem. (solution
-	      // component i vanished, and atol or atol(i) = 0.)
-    case -5:  // repeated convergence failures (perhaps bad jacobian
-	      // supplied or wrong choice of mf or tolerances).
-    case -4:  // repeated error test failures (check all inputs).
-    case -3:  // illegal input detected (see printed message).
-    case -2:  // excess accuracy requested (tolerances too small).
-      integration_error = 1;
-      return ColumnVector ();
-      break;
+      switch (istate)
+	{
+	case -13: // Return requested in user-supplied function.
+	case -6:  // error weight became zero during problem. (solution
+	          // component i vanished, and atol or atol(i) = 0.)
+	case -5:  // repeated convergence failures (perhaps bad jacobian
+	          // supplied or wrong choice of mf or tolerances).
+	case -4:  // repeated error test failures (check all inputs).
+	case -3:  // illegal input detected (see printed message).
+	case -2:  // excess accuracy requested (tolerances too small).
+	  integration_error = 1;
+	  break;
 
-    case -1:  // excess work done on this call (perhaps wrong mf).
-      working_too_hard++;
-      if (working_too_hard > 20)
-	{
+	case -1:  // excess work done on this call (perhaps wrong mf).
+	  working_too_hard++;
+	  if (working_too_hard > 20)
+	    {
+	      (*current_liboctave_error_handler)
+		("giving up after more than %d steps attempted in lsode",
+		 iwork.elem (5) * 20);
+	      integration_error = 1;
+	    }
+	  else
+	    {
+	      istate = 2;
+	      goto again;
+	    }
+	  break;
+
+	case 2:  // lsode was successful
+	  retval = x;
+	  t = tout;
+	  break;
+	  
+	default: // Error?
 	  (*current_liboctave_error_handler)
-	    ("giving up after more than %d steps attempted in lsode",
-	     iwork[5] * 20);
-	  integration_error = 1;
-	  return ColumnVector ();
+	    ("unrecognized value of istate returned from lsode");
+	  break;
 	}
-      else
-	{
-	  istate = 2;
-	  goto again;
-	}
-      break;
-
-    case 2:  // lsode was successful
-      break;
-
-    default:
-      // Error?
-      break;
     }
 
-  t = tout;
-
-  return x;
+  return retval;
 }
 
 #if 0

--- a/liboctave/LSODE.h
+++ b/liboctave/LSODE.h
@@ -128,7 +128,7 @@
   
   LSODE (const ColumnVector& state, double time, const ODEFunc& f);
 
-  ~LSODE (void);
+  ~LSODE (void) { }
 
   void force_restart (void);
 
@@ -157,14 +157,15 @@
   int integration_error;
   int restart;
   int method_flag;
-  int *iwork;
-  double *rwork;
+  Array<int> iwork;
+  Array<double> rwork;
   int istate;
   int itol;
   int itask;
   int iopt;
   int liw;
   int lrw;
+  int working_too_hard;
 
   friend int lsode_f (int *neq, double *t, double *y, double *ydot);
 

--- a/liboctave/NPSOL.cc
+++ b/liboctave/NPSOL.cc
@@ -183,45 +183,40 @@
 
   // Constraint stuff.
 
+  double bigbnd = infinite_bound ();
+
   Matrix clin = lc.constraint_matrix ();
   double *pclin = clin.fortran_vec ();
 
-  Array<double> aclow (n+nclin+ncnln);
-  double *clow = aclow.fortran_vec ();
-
-  Array<double> acup (n+nclin+ncnln);
-  double *cup = acup.fortran_vec ();
+  ColumnVector aclow (n+nclin+ncnln);
+  ColumnVector acup (n+nclin+ncnln);
 
   if (bnds.size () > 0)
     {
-      for (int i = 0; i < n; i++)
-	{
-	  clow[i] = bnds.lower_bound (i);
-	  cup[i] = bnds.upper_bound (i);
-	}
+      aclow.insert (bnds.lower_bounds (), 0);
+      acup.insert (bnds.upper_bounds (), 0);
     }
   else
     {
-      double huge = 1.0e30;
-      for (int i = 0; i < n; i++)
-	{
-	  clow[i] = -huge;
-	  cup[i] = huge;
-	}
+      aclow.fill (-bigbnd, 0, n-1);
+      acup.fill (bigbnd, 0, n-1);
     }
 
-  for (int i = 0; i < nclin; i++)
+  if (nclin > 0)
     {
-      clow[i+n] = lc.lower_bound (i);
-      cup[i+n] = lc.upper_bound (i);
+      aclow.insert (lc.lower_bounds (), n);
+      acup.insert (lc.upper_bounds (), n);
     }
 
-  for (int i = 0; i < ncnln; i++)
+  if (ncnln > 0)
     {
-      clow[i+n+nclin] = nlc.lower_bound (i);
-      cup[i+n+nclin] = nlc.upper_bound (i);
+      aclow.insert (nlc.lower_bounds (), n+nclin);
+      acup.insert (nlc.upper_bounds (), n+nclin);
     }
 
+  double *clow = aclow.fortran_vec ();
+  double *cup = acup.fortran_vec ();
+
   Array<double> ac (ncnln);
   double *c = ac.fortran_vec ();
 

--- a/liboctave/QPSOL.cc
+++ b/liboctave/QPSOL.cc
@@ -87,13 +87,8 @@
 
   double bigbnd = infinite_bound ();
 
-  double *pa = 0;
-  Matrix clin;
-  if (nclin > 0)
-    {
-      clin = lc.constraint_matrix ();
-      pa = clin.fortran_vec ();
-    }
+  Matrix clin = lc.constraint_matrix ();
+  double *pa = clin.fortran_vec ();
 
   ColumnVector bl (n+nclin);
   ColumnVector bu (n+nclin);
@@ -105,8 +100,8 @@
     }
   else
     {
-      bl.fill (-bigbnd);
-      bu.fill (bigbnd);
+      bl.fill (-bigbnd, 0, n-1);
+      bu.fill (bigbnd, 0, n-1);
     }
 
   if (nclin > 0)