%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  This file solves the odefile BR.m and plots V and gating variables.
%  Plots currents IK, Ix, Is, and INa.
%  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all

global K Na X SS  nn t_old AH BH

nn=0;
t_old=0;

% s0 = [-79.5, 0.0005, 0.948, 0.999999, 0.0045, 0.999999, 0.056, 0.000000225];
% s0 = [-40, 0.0005, 0.948, 0.999999, 0.0045, 0.999999, 0.056, 0.000000225];

% Steady state intial condition:
% s0=[-79.5473971945913   0.000477040103889804 0.938198741554261 ...
%       0.93046862277159  0.00456209737137637  0.99995196944918  ...
%       0.558508702370444 2.24818892404663e-07];

% s0=[ -79.5    0.000477040103889804 0.938198741554261 ...
%       0.93046862277159  0.00456209737137637  0.99995196944918  ...
%       0.0558508702370444 2.24818892404663e-07];

s0=[ -84.52025813590147,
	0.01106002680373,
	0.98747810709297,
	0.97456548026162,
	0.00298864284111,
	0.99998108550269,
	0.01125776352785,
	0.00000017866136];

t_end = 200;

tspan = [1:1:t_end];

 [T,S] = ode15s('BR',tspan, s0, odeset('maxstep',1));  

figure(11);                     

 subplot(4,2,1),plot(T,S(:,1)),xlabel('time(ms)'), ylabel('V');      
 subplot(4,2,2),plot(T,S(:,2)),xlabel('time(ms)'), ylabel('m');      
 subplot(4,2,3),plot(T,S(:,3)),xlabel('time(ms)'), ylabel('h');      
 subplot(4,2,4),plot(T,S(:,4)),xlabel('time(ms)'), ylabel('j');    
 subplot(4,2,5),plot(T,S(:,5)),xlabel('time(ms)'), ylabel('d');      
 subplot(4,2,6),plot(T,S(:,6)),xlabel('time(ms)'), ylabel('f');      
 subplot(4,2,7),plot(T,S(:,7)),xlabel('time(ms)'), ylabel('x');    
 subplot(4,2,8),plot(T,S(:,8)),xlabel('time(ms)'), ylabel('Cai');  


figure(12);

subplot(2,1,1);
plot(T,K,'r',T,X,'b',T,SS,'m');
xlabel('time(ms)');
ylabel('IK,Ix,Is');
legend('IK','Ix','Is')
subplot(2,1,2);
plot(T,Na,'g');
xlabel('time(ms)');
ylabel('INa');
%V_init = S(:,1); m_init = S(:,2); h_init = S(:,3); n_init = S(:,4);
%save initial_400b V_init m_init h_init n_init
