Index: mttroot/mtt/bin/trans/make_ode2odes
==================================================================
--- mttroot/mtt/bin/trans/make_ode2odes
+++ mttroot/mtt/bin/trans/make_ode2odes
@@ -7,10 +7,13 @@
 ###############################################################
 ## Version control history
 ###############################################################
 ## $Id$
 ## $Log$
+## Revision 1.64  2001/08/08 02:15:00  geraint
+## Rationalisation of solver code, beginning with algebraic solvers.
+##
 ## Revision 1.63  2001/08/07 04:39:24  geraint
 ## Consolidated dassl and residual functions.
 ##
 ## Revision 1.62  2001/08/01 22:14:32  geraint
 ## Bug fix for dassl.
@@ -768,22 +771,25 @@
 
 #ifdef STANDALONE
 ColumnVector
 Fmtt_residual (const ColumnVector &X,
 	      const ColumnVector &DX,
-	      double t)
+	      double t,
+	      int &ires)
 {
 #else // !STANDALONE
 DEFUN_DLD (mtt_residual, args, , "mtt_residual")
 {
     static ColumnVector	X  (MTTNX+MTTNYZ);
     static ColumnVector	DX (MTTNX+MTTNYZ);
     static double      	t;
+    static int		&ires;
 
     X  = args(0).${vector_value} ();
     DX = args(1).${vector_value} ();
     t  = args(2).double_value ();
+    ires = static_cast<int>(args(3).double_value ());
 #endif // STANDALONE
 
     static ColumnVector residual (MTTNX+MTTNYZ);
     static ColumnVector U (MTTNU+MTTNYZ);
     static ColumnVector u (MTTNU);
@@ -828,11 +834,29 @@
 
     for (register int i = 0; i < MTTNX; i++)
       residual (i) = dx(i) - DX(i);
     
     if (MTTNYZ > 0)
-      residual.insert (yz,MTTNX);
+      {
+	residual.insert (yz,MTTNX);
+
+	// XXX:
+	// DASSL needs residual to be dependent on Ui and Uidot
+	// mtt_dassl always sets the initial Ui to zero, so
+	// Ui - h*Uidot should equal zero BUT, we don't know h
+	static double t_old;
+	double step;
+	if (t != t_old)
+	  {
+	    step = t - t_old;
+	    t = t_old;
+	  }
+	else
+	  step = t;	
+	for (register int i = MTTNX; i < MTTNX+MTTNYZ; i++)
+	  residual(i) += X(i) - DX(i)*step;
+      }
 
 #ifdef STANDALONE
     return residual;
 #else // !STANDALONE
     return octave_value (residual);

Index: mttroot/mtt/lib/cc/mtt_dassl.cc
==================================================================
--- mttroot/mtt/lib/cc/mtt_dassl.cc
+++ mttroot/mtt/lib/cc/mtt_dassl.cc
@@ -12,24 +12,25 @@
 #endif  // OCTAVE_DEV
 
 
 #ifdef STANDALONE
 extern ColumnVector
-Fmtt_residual (const ColumnVector &X, const ColumnVector &DX, double t);
+Fmtt_residual (const ColumnVector &X, const ColumnVector &DX, double t, int &ires);
 #endif // STANDALONE
 
 
 ColumnVector
-mtt_residual (const ColumnVector &X, const ColumnVector &DX, double t)
+mtt_residual (const ColumnVector &X, const ColumnVector &DX, double t, int &ires)
 {
 #ifdef STANDALONE
-  return Fmtt_residual (X, DX, t);
+  return Fmtt_residual (X, DX, t, ires);
 #else // !STANDALONE
   static octave_value_list args, f;
   args(0) = octave_value (X);
   args(1) = octave_value (DX);
   args(2) = octave_value (t);
+  args(3) = octave_value (ires);
   f = feval ("mtt_residual", args, 1);
   return f(0).VECTOR_VALUE ();
 #endif // STANDALONE
 }