#!/usr/bin/env octave

function ys = beuler ( f, fy, t, y_initial, h )
%
%  function ys = beuler ( f, fy, t, y_initial, h )
%
%  Uses the backward Euler method  
%  to estimate Y, the solution of an ODE, at equally spaced points in 
%  the range T(1) to T(2).  
%
%  Euler's forward method is used as a predictor,
%  and Newton's method is used to seek a solution of Euler's backward formula.
%
%  The name of the derivative function is F. 
%  The name of the partial derivative function is FY.
%
TOL = 0.00001;

N = length(y_initial);

x(1) = t;
nstep = ( t(2) - t(1) ) / h;
y(:,1) = y_initial;

for i = 1 : nstep

  x(i+1) = x(i) + h;
  yp = y(:,i) + h * feval ( f, x(i), y(:,i) );

  it = 0;
  resmax = TOL + 1.0;

  while ( resmax > TOL & it < 10 ) 
    r = yp - ( y(:,i) + h * feval ( f, x(i+1), yp ) );
    rprime = eye(N,N) - h * feval ( fy, x(i+1), yp );
    yp = yp - r / rprime;
    resmax = max ( r );
    it = it + 1;
  end

  ys(:,i+1) = yp;

end

## 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,1000];
y = euler (@Carelli05, t, y0, 0.025);


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