Index: mttroot/mtt/bin/trans/ode2odes_m
==================================================================
--- mttroot/mtt/bin/trans/ode2odes_m
+++ mttroot/mtt/bin/trans/ode2odes_m
@@ -13,10 +13,13 @@
 ###############################################################
 ## Version control history
 ###############################################################
 ## $Id$
 ## $Log$
+## Revision 1.13  1998/05/14 08:05:10  peterg
+## Put back under RCS
+##
 ## Revision 1.12  1998/02/25 18:02:39  peterg
 ## Removed the argument passing stuff .
 ## Replaced by the simpar.m method.
 ##
 ## Revision 1.11  1997/08/29 07:56:54  peterg
@@ -83,13 +86,13 @@
  
   %Read in state
   x = $1_state;
 
   %Set the initial output
-  if ny>0
-    y = $1_odeo(x,0);
-  end; 
+  %if ny>0
+  %  y = $1_odeo(x,0);
+  %end; 
  
   %Read in simulation parameters
   $1_simpar;
   T = [0:DT:LAST];
 
@@ -99,48 +102,57 @@
   %Defaults
   if exist('T')==0
     T=[0:1:100]
   end;
 
- % if exist('x0')==0
- %   % x0 = zeros(nx,1);
- %   x0 = x;
- %  end;
+  if exist('METHOD')==0
+    METHOD = 'Euler'
+  end;
+
+  if exist('x')==0
+    x = zeros(nx,1);
+   end;
 
   [n,m]=size(T);
   if m>n
     T=T';
   end;
 
-if nx>0
-%  x = lsode('$1_ode', x0, T);
-
-%Euler integration
-%  x = x0;
-  X=[]; Y=[];
-  dt = (T(2)-T(1))/STEPFACTOR;
-
-  for t=T'
-    X = [X; x'];
-    Y = [Y; y'];
-    ts = t;
-    for i=1:STEPFACTOR
-      dx = $1_ode(x,ts);
-      ts = ts + dt;
-      x = x + dx*dt;
-      if ny>0
-        y = $1_odeo(x,ts);
-      end;
-    end;
-  end;
-
+method = tolower(METHOD)
+
+if nx>0
+  if strcmp(method,'lsode')
+    X = lsode('$1_ode', x, T);
+  elseif strcmp(method,'euler')
+    %Euler integration
+    X=[];
+    dt = (T(2)-T(1))/STEPFACTOR;
+    for t=T'
+      X = [X; x'];
+      ts = t;
+      for i=1:STEPFACTOR
+        dx = $1_ode(x,ts);
+        ts = ts + dt;
+        x = x + dx*dt;
+      end;
+    end;
+  else
+    error('Method %s not available here', METHOD);
+    return;
+  end;
   write_matrix([T,X], '$1_odes');
 else
   X = zeros(size(T));
 end;
 
-if ny>0
+if ny>0 % compute y and print it
+  i = 0; Y=[];
+  for t=T'
+    i = i+1;
+    y = $1_odeo(X(i,:),t);
+    Y = [Y; y'];
+  end;
   write_matrix([T,Y], '$1_odeso');
 end;
 
 
 EOF