function [ V, taxis, n_spikes ] = LIF( I, C, R, V0, Vth, Vreset, Vrest, N )
%Simulate an LIF neuron.
%
%INPUTS:
%  I = input current.
%  C = capacitance.
%  R = resistance.
%  V0 = intial voltage.
%  Vth = threshold voltage.
%  Vreset = voltage reset value.
%  Vrest = resting voltage.
%  N = # steps to simulate.
%
%OUTPUT:
%  V = voltage of IF model.
%  taxis = time axis of simulation.
%  n_spikes = # spikes detected.

dt=0.01;                      %Set the timestep.

n_spikes = 0;                 %The # of spikes detected.

Vstar = I*R + Vrest;          %Define Vstar = target voltage.
tau = R*C;                    %Define tau = time constant

V = zeros(N,1);
V(1)=V0;                      %Set the initial condition.

for k=1:N-1                                     %For each time step,
    V(k+1) = V(k) + dt*(-(V(k)-Vstar)/tau);     %...evaluate V at the next time step
    if V(k+1) > Vth                             %When we cross threshold
        V(k+1) = Vreset;                        %...set V to reset
        n_spikes = n_spikes+1;                  %...and add 1 to spike counter.
    end
end

taxis = (1:length(V))*dt;         %Return the time axis too.

end

