V_Cmin=min(x(:,1));V_Cmax=max(x(:,1));
I_Amin=min(x(:,2));I_Amax=max(x(:,2));
I_Kmin=min(x(:,4));I_Kmax=max(x(:,4));

V_Crange=abs(V_Cmax-V_Cmin);
I_Arange=abs(I_Amax-I_Amin);
I_Krange=abs(I_Kmax-I_Kmin);

%I_Amin=0;I_Amax=1.25;I_Arange=I_Amax-I_Amin;

range=max([V_Crange,I_Arange,I_Krange]);

rangeslack=.1;

lnwd=1.5;

%%%%%%%%%%%%%%%%%%%
% Plot V_C-nullsurface
%%%%%%%%%%%%%%%%%%%

clear a b V_C I_A I_K cl c2 c3

[a,b]=meshgrid(0:.025:1, 0:.025:1);
%[a,b]=meshgrid(0:.01:1, 0:.01:1);
c1=0.*a+0;c2=0.*a+0.25;c3=0.*a+0.6;

I_K=I_Kmin-rangeslack*I_Krange+a*(I_Kmax+2*rangeslack*I_Krange-I_Kmin);
V_C=V_Cmin-rangeslack*V_Crange+b*(V_Cmax+3*rangeslack*V_Crange-V_Cmin);
I_A=-(Na(V_C-E_Na)+I_K+Cl(V_C-E_Cl)-Iin);

mesh(I_A,V_C,I_K,c1);

hold on

%%%%%%%%%%%%%%%%%%%
% Plot I_K-nullsurface
%%%%%%%%%%%%%%%%%%%

I_Amin=min(min(I_A));I_Amax=max(max(I_A));
%I_Amin=-.39;I_Amax=-.32;
I_Arange=abs(I_Amax-I_Amin);

clear V_C I_A I_K 

I_A=I_Amin-rangeslack*I_Arange+a*(I_Amax+2*rangeslack*I_Arange-I_Amin);
I_K=I_Kmin-rangeslack*I_Krange+b*(I_Kmax+2*rangeslack*I_Krange-I_Kmin);
V_C= K(I_K)+E_K;

mesh(I_A,V_C,I_K,c2);
%xlabel('V_C'); ylabel('I_A'); zlabel('I_{Na}');
hold on;



%%%%%%%%%%%%%%%%%%%
% Plot I_A-nullsurface
%%%%%%%%%%%%%%%%%%%

clear V_C I_A I_K cl

I_A=I_Amin-rangeslack*I_Arange+a*(I_Amax+2*rangeslack*I_Arange-I_Amin);
I_K=I_Kmin-rangeslack*I_Krange+.5*b*(I_Kmax+2*rangeslack*I_Krange-I_Kmin);
V_C=R*I_A;

mesh(I_A,V_C,I_K,c3);

legend('{V_C}\prime = 0','I_K\prime = 0','{I_A}\prime = 0')
xlabel('I_A','fontsize',16)
ylabel('V_C','fontsize',16)
zlabel('I_K','fontsize',16)


axis([-.55 .22 -.5 .55 .0 0.57])

%%%%%%%%%%%%%%%%%%%
%%% Intersection of I_A-nullsurface with I_K-nullsurface
%%%%%%%%%%%%%%%%%%%
clear V_C I_A I_K 
%I_A = (I_Amin-2*rangeslack*I_Arange):0.1*I_Arange:(I_Amax+2*rangeslack*I_Arange);
I_A =I_Amin:0.1:I_Amax;
V_C=R*I_A;
I_K=(V_C-E_K)/r1;
plot3(I_A,V_C,I_K,'k-','linewidth',lnwd);


%%%%%%%%%%%%%%%%%%%
%%% Intersection of V_C-nullsurface with I_K-nullsurface
%%%%%%%%%%%%%%%%%%%
clear V_C I_A I_K 
%I_Kmin=-.4;I_Kmax=1.4
I_K=I_Kmin-.5*rangeslack*I_Krange:0.001:I_Kmax+.5*rangeslack*I_Krange;
V_C= K(I_K)+E_K;
I_A=-(Na(V_C-E_Na)+I_K+Cl(V_C-E_Cl)-Iin);
plot3(I_A,V_C,I_K,'b-','linewidth',lnwd);




%%%%%%%%%%%%%%%%%%%
% Drop-off, and Touch-up Line
%%%%%%%%%%%%%%%%%%%
clear V_C I_A I_K cl
I_A = (I_Amin-rangeslack*I_Arange):0.1*I_Arange:(I_Amax+rangeslack*I_Arange);
V_C = K(i2)+E_K;
I_K = (V_C-E_K)/r1;
%figure(1)
plot3(I_A,V_C+0*I_A,I_K+0*I_A,'m--', 'linewidth',lnwd);
clear V_C I_K
V_C = K(i1)+E_K;
I_K = (V_C-E_K-(-r1+r2)*(i2-i1))/r1;
plot3(I_A,V_C+0*I_A,I_K+0*I_A,'m--', 'linewidth',lnwd);
