function [ output_args ] = SteadyState( input_args )

%This subfunction is used to calculate the steady state values
%of the state variables and uses ode45 to integrate the equations.

global IStim GK CM EK GNa GL ENa EL WholeTime TTotal V INa IK IL;

%Model constants defined in main program

%Initial conditions
V0=-60;
%Initial Mss
AlphaM=(V0+40)/(10*(1-exp(-(V0+40)/10)));
BetaM=4*exp(-(V0+65)/18);
Mss0=AlphaM/(AlphaM+BetaM);
    
%Initial Hss
AlphaH=0.07*exp(-(V0+65)/20);
BetaH=1/(1+exp(-(V0+35)/10));
Hss0=AlphaH/(AlphaH+BetaH);
    
%Initial Nss
AlphaN=(V0+55)/(100*(1-exp(-(V0+55)/10)));
BetaN=0.125*exp(-(V0+65)/80);
Nss0=AlphaN/(AlphaN+BetaN);
    
%Use ode45 to integrate system.
%File dv.m specifies the derivatives of the state variables when
%a stimulus is applied.
%File dv0.m specifies these derivatives in the absence of a stimulus.
%T1 specifies the time during which stimulus is applied, T2 is
%post-stimulus simulation time.

  [T1,V1]=ode45('dv',[0,0.05],[-60; Mss0; Hss0; Nss0]);
[T2,V2]=ode45('dv0',[0.05,WholeTime],[V1(length(V1),1); V1(length(V1),2); V1(length(V1),3); V1(length(V1),4)]);
   
    
% Concatenate results with and without stimulus to get full result
TTotal=[T1;T2];
V=[V1(:,1);V2(:,1)];
M=[V1(:,2);V2(:,2)];
H=[V1(:,3);V2(:,3)];
N=[V1(:,4);V2(:,4)];
    
    
% Use state variable results to calculate membrane currents    
INa1=GNa*V1(:,2).^3.*V1(:,3).*(V1(:,1)-ENa);
IK1=GK*V1(:,4).^4.*(V1(:,1)-EK);
IL1=GL*(V1(:,1)-EL);
INa2=GNa*V2(:,2).^3.*V2(:,3).*(V2(:,1)-ENa);
IK2=GK*V2(:,4).^4.*(V2(:,1)-EK);
IL2=GL*(V2(:,1)-EL);
INa=[INa1;INa2];
IK=[IK1;IK2];
IL=[IL1;IL2];
    
ITotal1=INa1+IK1+IL1-IStim;
ITotal2=INa2+IK2+IL2;
ITotal=[ITotal1;ITotal2];
    
    
%Plot results: state variables (fig.1) and currents (fig.2)
     %    figure (1)
%V
subplot(1,1,1)
     plot (TTotal,V,'k')
     axis([0 WholeTime -80 60])
     ylabel('Membrane Voltage (mV)')
     xlabel('Time (msec)')

%M
% subplot(3,1,2)
     %plot (TTotal,M,'black')
     %hold on
%ylabel('M, H, N')

%H
%    plot (TTotal,H,'blue')
     %    text(0.8,0.75,'m')
     %    text(5,0.75,'n')
     %    text(7,0.2,'h')
    
%N
%plot (TTotal,N,'red')
%    hold off


%ITotal
% subplot(3,1,3)
     %plot (TTotal,ITotal,'black')
    
     %xlabel('Time')
     %ylabel('I_T_o_t_a_l')

% individual ionic currents
%    figure (2)
     %plot (TTotal,INa,'r') %INa - T
%hold on
%grid on
%plot (TTotal,IK,'black') %IK - T
%plot (TTotal,IL,'g') %IL - T
%title('I_N_a(red),I_K(black),I_L(green) - T')
     %xlabel('Time')
     %ylabel('I_N_a,I_K,I_L')
