function be_vs_trap

VCl = -68;      % mV
A = 4*pi*1e-6;  % cm^2 patch area
Cm = 1;         % micro F/cm^2
gCl = 0.3;      % mS/cm^2
tau = Cm/gCl;

dt = [.5 .1 .05 .01 .005 .001];

for k=1:6,

    v = bepsa(dt(k),20);

    cnv = trapsa(dt(k),20);

    t = 0:dt(k):20;

    vex = VCl + (2e-5)*exp(-t/tau).*(t.^2)/2/A/Cm/tau;

    %err(k) = max(abs(v-vex))/max(abs(vex));

    be_err(k) = max( abs(v-vex));

    cn_err(k) = max( abs(cnv-vex) );

end

figure(2)
loglog(dt,be_err,'o--','linewidth',2)
hold on
loglog(dt,cn_err,'x-','linewidth',2)
hold off
xlabel('dt  (ms)','fontsize',16)
ylabel('maximum absolute error  (mV)','fontsize',16)
legend('Backward Euler','Trapezoid','location','best')

%
% Backward Euler on a Passive Sphere
% with alpha-like input
%
%  usage   v = bepsa(dt,Tfin)
%
%  e.g.    v = bepsa(0.01,40)
%

function v = bepsa(dt,Tfin)


VCl = -68;      % mV
A = 4*pi*1e-6;  % cm^2 patch area
Cm = 1;         % micro F/cm^2
gCl = 0.3;      % mS/cm^2
tau = Cm/gCl;   % ms

j = 1;
v = zeros(1,Tfin/dt);
v(1) = VCl;
t(1) = 0;

while t < Tfin

   t = j*dt;

   Istim = (2e-5)*t*exp(-t/tau)/tau;

   v(j+1) = (v(j) + dt*(VCl/tau + Istim/A/Cm))/(1+dt/tau);

   j = j + 1;

end

%
% Trapezoid on a Passive Sphere with alpha-like input
%
%  usage   v = trapsa(dt,Tfin)
%
%  e.g.    v = trapsa(0.01,40)
%

function v = trapsa(dt,Tfin)


VCl = -68;      % mV
A = 4*pi*1e-6;  % cm^2 patch area
Cm = 1;         % micro F/cm^2
gCl = 0.3;      % mS/cm^2
tau = Cm/gCl;   % ms

j = 1;
v = zeros(1,Tfin/dt);
t = v;
v(1) = VCl;
Istim(1) = 0;

while t < Tfin

   t = j*dt;

   Istim(j+1) = (2e-5)*t*exp(-t/tau)/tau;

   v(j+1) = ( (1-dt/tau/2)*v(j) + dt*( VCl/tau + (Istim(j)+Istim(j+1))/2/A/Cm ) )/(1+dt/tau/2);

   j = j + 1;

end