function [y,t]=beuler(f,fy,t,y_initial,h)

## backward euler algorithm.
## inputs: N, number of steps
##         a,b: t interval in the form t=[a,b]
##         w, initial data in vector form (one row)
##         x: initial guess in vector form
## outputs: out: [x1,x2,...xn,t] all different values in columns

  TOL = 0.00001;
  N = length(y_initial);

  x = [t(1)];
  nstep = ( t(2) - t(1) ) / h;
  y(:,1) = y_initial;
  
  ## algorithm
  for i=1:N
    y(:,i+1)=newton(f,fy,x,y(:,i),h,TOL); ## solve for y(i+1)
    x(i+1)=x(i)+h; ## increase t at t+h
  endfor

endfunction

function y=newton(f,fy,x,yi,c,TOL)
## solve for n non linear set of equations
## calls f.m and jaco.m
## inputs: x initial guesses
##         wj: to be used for backward euler: w(j,:)
##         c: to be used for backward euler: h
## outputs: out=answer vector.in the form [x1,x2,...xn]

  it = 0;
  resmax = TOL + 1.0;
  while ( resmax > TOL & it < 10 ) 
    x
    yi
    fy(x,yi)
    y=-inv(fy(x,yi))*f(x,yi); ## solves jaco(x)*y=-f(x)
    x=y'+x; ## newton iteration (update value)
    resmax=(sum((y.^2)))^.5; ## calculates ||y||
    it = it+1
  endwhile
  
endfunction

## Hodgkin-Huxley model driver for Octave

##hh_defs = "hh-substrate.m";
hh_defs = "hh-superstrate.m";

autoload ("hodgkin_huxley", hh_defs );
autoload ("hodgkin_huxley_ddy", hh_defs );
autoload ("hodgkin_huxley_init", hh_defs );

hodgkin_huxley_init; 

y0 = hodgkin_huxley_init(-65);

#t = linspace(0,100,1000);
#lsode_options("initial step size",0.001);
#y = lsode ("hodgkin_huxley", y0, t);

#t = [0,100];
#P = odeset ('InitialStep', 0.0001 );
#[t, y] = odesx (@hodgkin_huxley, t, y0, P);

t = [0,10];
y = beuler (@hodgkin_huxley, @hodgkin_huxley_ddy, t, y0, 0.0005);

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