function V=run(eNa,eK,gNa,gK)
% usage:
% V0=run(50,-85,4,0.35);
% V1=run(50,-85,8,0.35); 
% t = 0:0.02:15:
% plot(t,V0,t,V1)
% plot(t(2:end),diff(V1)./diff(t),t(2:end),diff(V0)./diff(t))

if (nargin<4)
  disp('Not all constants given, will use defalt values:');
  eNa = 50	
  eK = -85
  gNa = 4
  gK = 0.35
end

T = 50; % total simulation time

n = 60; % number of grid points in each direction
m = 8; % number of state variables 
S = zeros(n,n,m); % the solution tensor, S(:,:,1) is the membrane potential
x = linspace(0,1,n);

% set up the initial value
y0 = [-84,0.5e-7,0,0,1,1,0,1];
for i = 1:m
   S(:,:,i) = y0(i)*ones(n,n);
end
Q = round(3*n/8):round(5*n/8); 
q = length(Q)
S(Q,Q,1) = -20*ones(q,q); %adds a peak in the membrane potenial in the center
mesh(S(:,:,1))

% Q = round(n/2):n; 
% q = length(Q)
% S(Q,Q,1) = -20*ones(q,q); %adds a peak in v in the upper right corner
% contourf(x,x,S(:,:,1)); caxis([-80,40]); colorbar

% R = 1:round(n/2); 
% r = length(R)
% S(R,:,6) = zeros(r,n); %sets gate j to zero for x<0.5
% figure(3)
% mesh(S(:,:,6))

figure(1)

% discretization paramters
dt = 0.02; % the time step
dx = 1/n; % the space step (the domain is [0,1]^2)
d = 0.001; % the conductivity
alpha = d*dt/dx^2;

V = []; % The transmembrane potential in the corner of the domain
i = 0;
plot_dt = 1.0 % will only plot the solution for every plot_dt
t = 0:dt:T;
for i = 1:length(t)
   S = S + dt*ds_dt(t(i), S,[],eNa,eK,gNa,gK); % Solve the ODEs
   S(:,:,1) = heat(S(:,:,1),alpha);	% Solves the PDE
   V(i) = S(1,1,1); % keep the corner value
   if rem(t(i),plot_dt)<dt
     t(i)
     contourf(x,x,S(:,:,1)); caxis([-80,40]); colorbar
     %mesh(x,x,S(:,:,1)); axis([0 1 0 1 -90 30])
     pause(0.01)
   end
	
end
figure(2)
plot(t,V);  xlabel('ms'); ylabel('mV');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function V = heat(V, alpha)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[n,m] = size(V);

% Neumann condtions on the boundary
I = 1:n;
Ip = [2:n, n-1];
Im = [2,1:n-1];
J = 1:m;
Jp = [2:m, m-1];
Jm = [2,1:m-1];

% Explicit update
V(I,J) = V(I,J) + alpha*(V(Ip,J) + V(Im,J) - 4*V(I,J) + V(I,Jp) + V(I,Jm));


