Why do i receive error message "index exceeds number of array elements" here?

clear
vthresh = -50;
vreset = -65;
tmax = 2;
vleak = -70;
rmem = 100;
cmem = 0.1;
delta_t = 0.0001;
iappvec = 100:100:600;
tvec = 0:delta_t:tmax;
ref_period = 2.5;
for x = 1:length(iappvec) %cycling through iapp values
    vmem=zeros(1,length(tvec));
    vmem(1)=vreset;
    for i = 2:length(tvec) %forward euler
        dqdt=(vleak+vmem(i-1))/rmem+iappvec(x);
        dvdt=dqdt/cmem;
        vmem(i)=vmem(i-1)+dvdt*delta_t;
        spiketime=0;
        if vmem>vthresh %spike
            vmem=vreset;
            spiketime=tvec(i);
        end
        if tvec(i)<spiketime+ref_period
            vmem=vreset;
        end
    end
end

Index exceeds the number of array elements. Index must not exceed 1.

Answer 1

[Image Analyst]

vmem is a single number, not a vector

vmem(1)=vreset; % vmem = -65

So when you do

dqdt=(vleak+vmem(i-1))/rmem+iappvec(x);

and try to reference vmem(3-1) = vmem(2), well, there is no vmem(2).

Just step through it with the debugger to prove it.

Debugging in MATLAB | Doug's MATLAB Video Tutorials

[Walter Roberson]

The code had

    vmem=zeros(1,length(tvec));

so vmem was not a scalar at the time of

    vmem(1)=vreset; % vmem = -65

However, inside if statements, the code overwrites all of vmem, turning it into a scalar.

[Image Analyst]

Yes, it was a scalar on the first iteration (which was not a problem then), then it made it into a vector. Then the if statement made it back into a scalar again so that when the third iteration rolls around it's back to a scalar, thus causing the error.

Maybe a fix would be

vmem = vreset * ones(1,length(tvec)); % Still a vector after this, not a scalar.

[Walter Roberson]

? Where is vmem a scalar before it is assigned a vector by

vmem=zeros(1,length(tvec));

?

The fix is very likely to change

vmem=vreset;

to

vmem(i)=vreset;

[Image Analyst]

You're right. The first time through it was a vector with 20,001 elements, all zero except for the first one which was -65. The second time through, it was a scalar with a single element/value of -65. The lack of comments makes it hard to figure out what's going on, like when/if you need to reset vmem.

Answer 2

[Walter Roberson]

for x = 1:length(iappvec) %cycling through iapp values
    vmem=zeros(1,length(tvec));

vmem is a vector reinitialized for each x value.

        if vmem>vthresh %spike
            vmem=vreset;
            spiketime=tvec(i);
        end
        if tvec(i)<spiketime+ref_period
            vmem=vreset;
        end

In both of those conditions, you overwrite all of the vector vmem, making it into a scalar.