function dydt = axon_deriv ( t, y )

%*****************************************************************************80
%
%% axon_deriv() evaluates the right hand side of the axon ODE.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    22 April 2021
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Alan Hodgkin, Andrew Huxley,
%    A quantitative description of membrane current and its
%    application to conduction and excitation in nerve,
%    The Journal of Physiology,
%    Volume 117, Number 4, pages 500-544, August 1952.
%
%  Input:
%
%    real T, Y(4): the time and solution value.
%
%  Output:
%
%    real DYDT(4): the derivative value.
%

%
%  Extract individual variable values at time t.
%
  V = y(1);
  n = y(2);
  m = y(3);
  h = y(4);
%
%  Get parameter values.
%
  [ C, E_K, E_Na, E_rest, G_K, G_Na, t0, y0, tstop ] = axon_parameters ( );
%
%  Compute alpha and beta for n:
%  Equations 12 and 13.
%
  alpha_n = 0.01 * ( ( 10.0 - V ) / ( exp ( ( 10.0 - V ) / 10.0 ) - 1.0 ) );
  beta_n = 0.125 * exp ( - V / 80.0 );
%
%  Compute alpha and beta for m:
%  Equations 20 and 21.
%
  alpha_m = 0.1 * ( ( 25.0 - V ) / ( exp ( ( 25.0 - V ) / 10.0 ) - 1.0 ) );
  beta_m = 4.0 * exp ( - V / 18.0 );
%
%  Compute alpha and beta for h:
%  Equations 23 and 24.
%
  alpha_h = 0.07 * exp ( - V / 20.0 );
  beta_h = 1.0 / ( exp ( ( 30.0 - V ) / 10.0 ) + 1.0 );
%
%  Calculate the currents.
%  I_Na:  equations 3 and 14.
%  I_K:   equations 4 and 6.
%  I_ion: equation 26.
%
  I_Na = m^3 * G_Na * h * ( V - E_Na );
  I_K = n^4 * G_K * ( V - E_K );
  I_ion = axon_current ( t ) - I_K - I_Na;
%
%  Evaluate the derivatives.
%
%  dndt: equation 7
%  dmdt: equation 15
%  dhdt: equation 16
%
  dVdt = I_ion / C;
  dndt = alpha_n * ( 1.0 - n ) - beta_n * n;
  dmdt = alpha_m * ( 1.0 - m ) - beta_m * m;
  dhdt = alpha_h * ( 1.0 - h ) - beta_h * h;
%
%  Pack derivatives into vector.
%
  dydt = [ ...
    dVdt; ...
    dndt; ...
    dmdt; ...
    dhdt ];

  return
end

