function out = neuron( task, opt)

% Declaration of global variables to be used repeatedly
% any variable not declared as "global" will be cleared after each call
global    model               % 
global	ts ns vname		% time span, state dim, and variable names
global	Inp				% current input
global	T Xt It 		      % time, state, and input tables
global    gna gk ga gca gahp  % conductance
% If you want to check or change these variables from command line, type, e.g.
% >>global ts

% If you just type
% >>neuron
% print Usage
if nargin == 0
   error('Usage : neuron [init/run]')
   return
end
% If a task is specified:
switch(task)
case 'init'				% Initialization
   if nargin >= 2, model=opt, end;
   if isempty(model), model='cclamp', end;
   figure(1);clf;
   % default time span and initial state
   [ts,x] = feval(model,[],[],'init');
   % is equivalent to, e.g., cclamp([],[],'init');
   ns = size(x,1);		% state space dimension
   vname = feval( model,[],[],'name');	% variable names
   neuron( 'input',10);	% default stimulus
   T = 0;				% initialize time, state, input trace, and conductance
   Xt = x';
   It = 0;
   gna = 221;
   gk = 47;
   ga = 0;
   gca = 8.5;
   gahp = 7;
case 'run'				% Run the simulation
   % options for ode simulation:
   odeopt = odeset('OutputFcn','');
   % integration of 'stiff' differential equation
   % time steps are chosen automatically
   % Inp could be constant or (t,I) table
   if size(Inp,2)==2, ts=Inp([1,end],1); end
   [ T1, X1] = ode23s(model, ts, Xt(end,:)',[],Inp,gna,gk,ga,gca,gahp);
   % concatenate with previous output
   T = [ T; T(end)+T1];	% time 
   Xt = [ Xt; X1];		% state
   if size(Inp,2)==2	% input
      It = [ It; interp1( Inp(:,1), Inp(:,2), T1)];	% table
   else
      It = [ It; repmat( Inp(1), size(T1))];		% const
   end
   figure(1); neuron('wave');
   
case 'wave'			% Plot waveforms
   for i = 1:ns
      subplot(ns+1,1,i);	% plot i-th state variable
      plot( T, Xt(:,i));
     ylabel( vname(i,:));
   end
   subplot(ns+1,1,ns+1);	% plot current input
   plot( T, It(:));
   ylabel( 'I'); xlabel( 't (ms)');
   
case 'input'              % Input
   if ischar(opt)
      opt = eval(opt);
   else
      opt = 10;
   end	% strig to value
   Inp = [0 0;20 0;20.1 opt;100 opt]
   
case 'save'			% Save all global variables, e.g.: >>cm save 'data1'
   save( opt);
   
case 'load'			% Relaod all global variables, including T, X
   load( opt);
   
otherwise
   error( [' invalid task: ', task]);
   
end