%2014-01-23

% Passive membrane equation

clear all
close all

C = 1;
El = -65;
Gl = 0.1;
Iapp = 1;

T = 1000;
dt = 0.1;
t = 0:dt:T;

V = zeros(1,length(t));

ti = 100;
tf = 400;
H = zeros(1,length(t));
H(floor(ti/dt):floor(tf/dt))=1;

V(1) = -65;

for j=1:length(t)-1
    kv1 = (-Gl*(V(j)-El)+Iapp*H(j))/C;
    av = V(j)+kv1*dt;
    kv2 = (-Gl*(av-El)+Iapp*H(j+1))/C;
    V(j+1) = V(j) + (kv1+kv2)*dt/2;
end

figure
hold on
plot(t,V,'b','linewidth',2);
axis([0 T -70 -40]);
set(gca,'fontsize',20);
xlabel('t')
ylabel('V');