Please help. Error of index out of bounds

function Y = Linda

N=100; % Total size

en=50; % plot every nth time interval

T=zeros(N+1,N+1); % T is the transition matrix, defined below

v=linspace(0,N,N+1);

t(1)=1;

t=1;

beta=0.23;

v=20;

gamma=1;

b=11;

p=zeros(t(1)+1,N+1);

%p(1,3)=1; % Two individuals initially infected.

bt=beta*v.*(N-v)/N; dbstop if error

dt=(b+gamma)*v;

for i=1:N % Define the transition matrix

T(i,i)=1-bt(i)-dt(i); % diagonal entries

T(i,i+1)=dt(i+1); % superdiagonal entries

T(i+1,i)=bt(i); % subdiagonal entries

end

T(1,1)=1;

T(1,2)=dt(2);

T(N+1,N+1)=1-dt(N+1);

for t=1:t(1)

y=T*p(t,:);

p(t+1,:)=y;

end

pm(1,:)=p(1,:);

for t=1:t(1)/en;

pm(t+1,:)=p(en*t,:);

end

ti=linspace(0,t(1),t(1)/en+1);

st=linspace(0,N,N+1);

mesh(st,ti,pm);

xlabel('Number of Infectives');

ylabel('Time Steps');

zlabel('Probability');

view(140,30);

axis([0,N,0,time,0,1]);

%%%%%%%%%%%%%%%%%%%%%

%This is the error i'm having--

%Attempted to access dt(2); index out of bounds because numel(dt)=1.

Error in Linda (line 18)

T(i,i+1)=dt(i+1); % superdiagonal entries

[Jan]

Please format your code properly to improve the readability. Then explain, what the problem is. "With some errors" is not clear. It is much easier to fix a problem than to guess, what you consider as a problem.

Answer 1

[Manvi Goel]

The following errors were there in your code:

1. inappropriate spaces in for loops
2. spaces between clear and the variable associated

I have fixed them and the code below now runs fine.

function y = Stoch
% A program for iterations
% y1 y2 y3 y4 and y5 are the different classes
% y10 y20 y30 y40 and y50 are the initial conditions
% Problem-dependent statements are marked with a %***
% icase=1 corresponds to the deterministic problem
% nt is the number of steps
% h is the step size
% Accuracy generally increases as h decreases
clear
for icase = 1:2
    clear tt
    clear yp1
    clear yp2
    clear yp3
    clear yp4
    clear yp5
    nsamp=100; %***
    tmax=50; %***
    nt=500; %***
    y10=20; %***
    y20=25; %***
    y30=15; %***
    y40=25; %***
    y50=15; %***
    
    if(icase==1) 
        nsamp=1; 
    end
    h=tmax/nt;
    hs=sqrt(h);
    randn('state',20); %initiates the random number generator
    te1=zeros(nsamp,1);
    te2=zeros(nsamp,1);
    te3=zeros(nsamp,1);
    te4=zeros(nsamp,1);
    te5=zeros(nsamp,1);
    te6=zeros(nsamp,1);
    jj1=0;
    jj2=0;
    jj3=0;
    jj4=0;
    jj5=0;
    jj6=0;
    for jj=1:nsamp
        y1=y10;
        y2=y20;
        y3=y30;
        y4=y40;
        y5=y50;
        yp1(1)=y1;
        yp2(1)=y2;
        yp3(1)=y3;
        yp4(1)=y4;
        yp5(1)=y5;
        r=randn(nt+1,14);
        nchk1=0;
        nchk2=0;
        nchk3=0;
        nchk4=0;
        nchk5=0;
        n=0;
        t=0;
        chk=0;
        tt(1)=0;
        while (chk==0)
            n=n+1;
            t=t+h;
            if(jj==nsamp) 
                tt(n+1)=t; 
            end
            Lambda=0.070;
            mu=0.0048;
            alpha=0.03;
            beta=0.01;
            delta=0.2;
            q_1=0.7;
            q_2=0.63;
            gamma_1=0.17;
            gamma_2=0.2;
            lambda_1=0.71;
            lambda_2=0.82;
            k=0.5;
            f1=Lambda-beta*y3*y1-mu*y1;
            f2=beta*y3*y1-(k+mu)*y2;
            f3=(1-delta)*k*y2-(alpha+gamma_1+mu)*y3;
            f4=delta*k*y2+alpha*y3-(gamma_2+mu)*y4;
            f5=(1-q_1)*gamma_1*y3+(1-q_2)*gamma_2*y4-mu*y5;
            g1=sqrt(delta)*r(n,1)-sqrt(mu*y1)*r(n,2)-sqrt(beta*y1*y3)*r(n,1);
            g2=sqrt(beta*y1*y3)*r(n,4)-sqrt(mu*y2)*r(n,5)-sqrt((1-delta)*k)*r(n,6)-sqrt(delta*k)*r(n,7); dbstop if error
            g3=sqrt((1-delta)*k)*r(n,5)-sqrt(mu*y3)*r(n,7)-sqrt(alpha*y3)*r(n,8)-sqrt((1-q_1)*gamma_1)*r(n,9);
            g4=sqrt(delta*k)*r(n,6)+sqrt(alpha*y3)*r(n,8)-sqrt(mu*y4)*r(n,11)-sqrt((1-q_2)*lambda_2)*r(n,12)-sqrt(q_2*gamma_2)*r(n,13);
            g5=sqrt((1-q_1)*gamma_1)*r(n,9)+sqrt(2-q_2*gamma_2)*r(n,13)-sqrt(mu*y5)*r(n,14);

            if(icase==1) 
                g1=0; 
            end
            if(icase==1) 
                g2=0; 
            end
            if(icase==1) 
                g3=0; 
            end
            if(icase==1) 
                g4=0; 
            end
            if(icase==1) 
                g5=0; 
            end

            y1=y1+h*f1+hs*g1;
            y2=y2+h*f2+hs*g2;
            y3=y3+h*f3+hs*g3;
            y4=y4+h*f4+hs*g4;
            y5=y5+h*f5+hs*g5;

            if(jj==nsamp) 
                yp1(n+1)=y1; 
            end
            if(jj==nsamp) 
                yp2(n+1)=y2; 
            end
            if(jj==nsamp) 
                yp3(n+1)=y3; 
            end
            if(jj==nsamp) 
                yp4(n+1)=y4; 
            end
            if(jj==nsamp) 
                yp5(n+1)=y5; 
            end
% This is Euler's approximation to the SDE
            if (y1 < 1)
                chk=1;
                jj1=jj1+1;
                te1(jj1)=t;
            end
            if (y2 < 1)
                chk=1;
                jj2=jj2+1;
                te2(jj2)=t;
            end
            if (y3 < 1)
                chk=1;
                jj3=jj3+1;
                te3(jj3)=t;
            end
            if (y4 < 1)
                chk=1;
                jj4=jj4+1;
                te4(jj4)=t;
            end
            if (y5 < 1)
                chk=1;
                jj5=jj5+1;
                te5(jj5)=t;
            end
            if (t >tmax)
                chk=1;
                jj6=jj6+1;
                te6(jj6)=t;
                chk=1;
            end
        end% end of while (chk==0) loop
    end% end of for jj=1:nsamp loop
    tp=0; tp1=0; tp2=0; tp3=0; tp4=0; tp5=0; tp6=0;
    if(jj1 ~= 0) 
        tp1=sum(te1)/jj1; 
    end
    if(jj2 ~= 0)
        tp2=sum(te2)/jj2;
    end
    if(jj3 ~=0) 
        tp3=sum(te3)/jj3;
    end
    if(jj4 ~= 0) 
        tp4=sum(te4)/jj4; 
    end
    if(jj5 ~= 0) 
        tp5=sum(te5)/jj5; 
    end
    if(jj6 ~= 0) 
        tp6=sum(te6)/jj6; 
    end

    if(jj1+jj2+jj3+jj4+jj5~=0)
        tp=(sum(te1)+sum(te2)+sum(te3)+sum(te4)+sum(te4))/(jj1+jj2+jj3+jj4+jj5); 
    end

    p1=jj1/nsamp;
    p2=jj2/nsamp;
    p3=jj3/nsamp;
    p4=jj4/nsamp;
    p5=jj5/nsamp;
    p6=jj6/nsamp;
    disp(' ')

    if(icase==1) 
        disp(' Deterministic Calculational Results'); 
    end
    if(icase==2) 
        disp(' Stochastic Calculation Results');
    end

    disp(' icasensamp h tmax')
    disp((sprintf(' %12.0f %12.0f %12.5f %12.2f',icase,nsamp,h,tmax)));
    disp(' tp1 p1')
    disp((sprintf(' %12.6f %12.6f', tp1, p1)));
    disp(' tp2 p2')
    disp((sprintf(' %12.6f %12.6f', tp2, p2)));
    disp(' tp3 p3')
    disp((sprintf(' %12.6f %12.6f', tp3, p3)));
    disp(' tp4 p4')
    disp((sprintf(' %12.6f %12.6f', tp4, p4)));
    disp(' tp5 p5')
    disp((sprintf(' %12.6f %12.6f', tp5, p5)));
    disp(' tp6 p6')
    disp((sprintf(' %12.6f %12.6f', tp6, p6)));
    disp(' tp p1+p2+p3+p4+p5')
    disp((sprintf(' %12.6f %12.6f', tp, p1+p2+p3+p4+p5)));

    if(icase==1) 
        title('Deterministic'); 
    end
    if(icase==2) 
        title('Stochastic'); 
    end

    set(gca,'fontsize',18,'linewidth',1.5);
    plot(tt,yp1,'r-')%,tt,yp2,'k-', tt,yp4,'r-',tt,yp3,'y-',tt,yp5,'g-')
    xlabel('Time t')
    ylabel(' POPULATIONS')
    legend('Infected'),% 'Exp','Inf','Qua','Rec')

    if(icase==2) 
        title('Stochastic'); 
    end
    hold on
end% end of for icase=1:2 loop
    hold off