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;