function varargout=hhode(t, y, flag)% HHODE - ODE file for HH equations. State vector is [V,n,m,h]'. External% current injected by function iexthh(t). Model parameters returned using% hhpars = [gnamx; Ena; gkmx; Ek; gleak; El; C] = gethhparms.% flag='': returns ydot (column)%     ='init': returns [tsp,y0,opt], time span, initial values, options.if nargin<3							% Catch for numjac the dumbjack	varargout{1} = hhdot(t,y);else	switch flag		case ''                % Return ydot vector			varargout{1} = hhdot(t,y);		case 'init'             % Return default tspan, y0, and options			[varargout{1:3}] = init;	%	case 'jacobian'         % Return Jacobian matrix	%		varargout{1} = hhjac(t,y);		otherwise			error(strcat('***Unknown flag=''',flag,''' in HHODE.***'))	endendreturnfunction ydot = hhdot(t,y)% Note the external current is provided by function iexthh(t)% Parameters (mS/cm^2, mV, uFd/cm^2)% hps = [gnamx=120; Ena=55; gkmx=36; Ek=-72; gleak=0.3; El=-49.4; C=1]    hps = gethhparms;	ydot=zeros(4,1);	% Membrane potential:	ydot(1) = (iexthh(t) - hps(1)*y(3)^3*y(4)*(y(1)-hps(2)) - ...	      hps(3)*y(2)^4*(y(1)-hps(4)) - hps(5)*(y(1)-hps(6)))/hps(7);% n	ydot(2) = (1-y(2))*alphn(y(1)) - y(2)*betan(y(1));% m	ydot(3) = (1-y(3))*alphm(y(1)) - y(3)*betam(y(1));% h	ydot(4) = (1-y(4))*alphh(y(1)) - y(4)*betah(y(1));returnfunction [tsp,y0,opt]=init   tsp = [0; 50];   y0 = [-60., 0.318, 0.053, 0.596]';   opt = odeset('OutputFcn','odeplot', 'jacobian','off', 'RelTol',1e-4);return