 clear all
 global iK iK1 iNa iSi ib iKp KKo1 Tm Th Tj Td Tf Tx step t_old nn t_end

 time1 = cputime;

 t_old = 0;
 nn = 0;
 step  = 2;
 t_end = 2000; %Brad had 8000 here, which is too big for my version of Matlab

%s0 = [-84.5527007229, 0.0016649036, 0.9833233322, 0.9895346980,...
%       0.0029760116, 0.9999812598, 0.0056403099, 0.0001783547];
%s0 = [40, 0.0016649036, 0.9833233322, 0.9895346980, 0.0029760116,...
%       0.9999812598, 0.0056403099, 0.0001783547];
%s0=[-84.55270072288010, 0.00166490359515, 0.98332333217090, 0.98953469798485,   0.00297601158336, 0.99998125984966, 0.00564030994985, 0.00017835466557,...
%    -84.55270072288010, 0.00166490359515, 0.98332333217090, 0.98953469798485,   0.00297601158336, 0.99998125984966, 0.00564030994985, 0.00017835466557];
s0 = [-30,0.00168,0.98334,0.9895,0.002977,0.99999,0.0056,0.0001785,...
       -30,0.00168,0.98334,0.9895,0.002977,0.99999,0.0056,0.0001785]; 
       %Ko = 5.4      
%s0=[-45.6511,0.4042,0.0065,0.0066,0.0966,0.9720,0.1726,0.0017],%Ko = 25.4

 tspan = [1:step:t_end+step];               
 [T,S] = ode15s('LR1_coupODE',tspan, s0, odeset('maxstep',1));  

%(S(t_end,:)-S(t_end-1,:));

figure(1);
subplot(2,1,2), plot(T,S(:,1),'b',T,S(:,9),'r','Linewidth',[2]), 
xlabel('t(ms)','FontSize',[14]);
ylabel('V(mV)','FontSize',[14]);
% subplot(3,2,1),plot(T,S(:,1)), ylabel('V');      
% subplot(3,2,2),plot(T,S(:,2)), ylabel('n');      
set(gca,'FontSize',[14]);
title('Coupled LR1 V vs. T for m=0.9, d = 0.1 and p=Ko = 30','FontSize',[14]);
legend('V1','V2',-1);

%figure;                     
% plot(T,S(:,1)),hold on,plot(T,S(:,9),'r'), ylabel('V');      

pp = 0;
if pp == 1
figure(2);                     

 subplot(4,2,1),plot(T,S(:,1)),hold on,plot(T,S(:,9),'r'), ylabel('V');      
 subplot(4,2,2),plot(T,S(:,2)),hold on,plot(T,S(:,10),'r'), ylabel('m');      
 subplot(4,2,3),plot(T,S(:,3)),hold on,plot(T,S(:,11),'r'), ylabel('h');      
 subplot(4,2,4),plot(T,S(:,4)),hold on,plot(T,S(:,12),'r'), ylabel('j');    
 subplot(4,2,5),plot(T,S(:,5)),hold on,plot(T,S(:,13),'r'), ylabel('d');      
 subplot(4,2,6),plot(T,S(:,6)),hold on,plot(T,S(:,14),'r'), ylabel('f');      
 subplot(4,2,7),plot(T,S(:,7)),hold on,plot(T,S(:,15),'r'), ylabel('x');    
 subplot(4,2,8),plot(T,S(:,8)),hold on,plot(T,S(:,16),'r'), ylabel('Cai');  


figure(3);                     
 subplot(3,2,1),plot(T,Tm),xlabel('time(ms)'), ylabel('Tau_m');      
 subplot(3,2,2),plot(T,Th),xlabel('time(ms)'), ylabel('Tau_h');      
 subplot(3,2,3),plot(T,Tj),xlabel('time(ms)'), ylabel('Tau_j');      
 subplot(3,2,4),plot(T,Td),xlabel('time(ms)'), ylabel('Tau_d');    
 subplot(3,2,5),plot(T,Tf),xlabel('time(ms)'), ylabel('Tau_f');      
 subplot(3,2,6),plot(T,Tx),xlabel('time(ms)'), ylabel('Tau_x');      
 figure;
 plot(T,KKo1),xlabel('time'),ylabel('External Potassium');
 figure;
 subplot(2,2,1),plot(T,iK,'r',T,iK1+iKp+ib,'k'),axis([0,t_end,-1,7]),...
 xlabel('time(ms)'),ylabel('iK,(iK1,iKp,ib)'),title('black = sum');
 subplot(2,2,1),plot(T,ib,'r',T,iK1+iKp+iK,'k'),axis([0,t_end,-1,7]),...
 xlabel('time(ms)'),ylabel('ib,(iK1,iKp,iK)'),title('black = sum');
 subplot(2,2,2),plot(T,iK,'r',T,iK1,'b',T,iKp,'m',T,ib,'k'),...
 axis([0,t_end,-1,7]),...
 xlabel('time(ms)'),ylabel('iK,iK1,iKp,ib'),...
 title('blue=iK1,red=iK,pur=iKp,black=ib');
 subplot(2,2,3),plot(T,iNa,'g'),xlabel('time(ms)'),ylabel('iNa');
 subplot(2,2,4),plot(T,iSi,'g'),xlabel('time(ms)'),ylabel('iSi');
 end
 time = cputime-time1
 %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
 
