#!/usr/bin/env octave

function ys = euler(f,t,y,h)
  t0=t(1); tn=t(2);
  tx=t0; ys=zeros(ceil((tn-t0)/h),length(y));
  ys(1,:)=y; i=1;
  while (tx<tn)
    i = i+1; tx=tx+h;
    ys(i,:) = ys(i-1,:)' + (h * feval(f,tx,ys(i-1,:)));
  endwhile
endfunction

## Carelli 05 model driver for Octave

carelli05_defs = "carelli05.m";

autoload ("Carelli05", carelli05_defs );
autoload ("Carelli05_init", carelli05_defs );

y0 = Carelli05_init(-68)

#t = linspace (0, 100, 1000);
#lsode_options("initial step size",1e-3);
#lsode_options("minimum step size",1e-5);
#lsode_options("absolute tolerance",1e-2);
#lsode_options("relative tolerance",1e-4);
#y = lsode ("Carelli05", y0, t);

t = [0,100];
y = euler (@Carelli05, t, y0, 1e-3);


save "-ascii" "carelli05.dat" y;