i need to find error in code because I it's increasing but it should be decresing

i need to find error in code because I it's increasing but it should be decresing

the equation of Sir model in runge kutta (hune's method)

this is laws

i get the constants on this picture and laws

this my code

h=1.5;
a=4*10^-6; %a=alfa
y=0.017; %gamaa
S=zeros(22,1);
S(1)=54000;
I=zeros(22,1);
I(1)=757;
R=zeros(22,1);
R(1)=0;
fprintf('\n')
disp('                                      heun`s method ')
disp('-----------------------------------------------------------------------------------------------')
disp('    t              S                 I                   R')
disp('-----------------------------------------------------------------------------------------------')
fprintf('\n')
for i=1:22
    k1=-a*S(i)*I(i);
    l1=a*S(i)*I(i)-y*I(i);
    m1=y*I(i);
    k2=-a*(S(i)+h*k1)*(I(i)+h);
    l2=a*(S(i)+h*k1)*(I(i)+h)-y*(I(i)+h);
    m2=y*(I(i)+h);
    S(i+1)=S(i)+(h/2)*(k1+k2);
    I(i+1)=I(i)+(h/2)*(l1+l2);
    R(i+1)=R(i)+(h/2)*(m1+m2);
    
    fprintf('%10f %17f %17f %17f\n',i,S(i+1),I(i+1),R(i+1))
end
t=[1:22];
hold on
plot(S(t),t,'-b')
plot(I(t),t,'--r')
plot(R(t),t,'-g')

Answer 1

[Mathieu NOE]

hi

seems you have swapped the x and y data when you do the plot

after the fix we get this :

full code a bit improved :

h=1.5;
alfa=4*10^-6; %a=alfa
gamma=0.017; %gamma
n = 22; % time steps (samples)
S=zeros(n,1);
I=zeros(n,1);
R=zeros(n,1);
S(1)=54000;
I(1)=757;
R(1)=0;
fprintf('\n')
disp('                                      heun`s method ')
disp('-----------------------------------------------------------------------------------------------')
disp('    t              S                 I                   R')
disp('-----------------------------------------------------------------------------------------------')
fprintf('\n')
for i=1:n
    k1=-alfa*S(i)*I(i);
    l1=alfa*S(i)*I(i)-gamma*I(i);
    m1=gamma*I(i);
    
    k2=-alfa*(S(i)+h*k1)*(I(i)+h);
    l2=alfa*(S(i)+h*k1)*(I(i)+h)-gamma*(I(i)+h);
    m2=gamma*(I(i)+h);
    
    S(i+1)=S(i)+(h/2)*(k1+k2);
    I(i+1)=I(i)+(h/2)*(l1+l2);
    R(i+1)=R(i)+(h/2)*(m1+m2);
    
    fprintf('%10f %17f %17f %17f\n',i,S(i+1),I(i+1),R(i+1))
end
t=[1:n];
plot(t,S(t),'-b',t,I(t),'--r',t,R(t),'-g')
legend('S','I','R');

Answer 2

[James Tursa]

You need to use all the derivatives at the first point to propagate for the initial guess at the second point. I.e., these lines

k2=-a*(S(i)+h*k1)*(I(i)+h);
l2=a*(S(i)+h*k1)*(I(i)+h)-y*(I(i)+h);
m2=y*(I(i)+h);

should be this instead, where (I(i)+h) is changed to (I(i)+h*l1):

k2=-a*(S(i)+h*k1)*(I(i)+h*l1);
l2=a*(S(i)+h*k1)*(I(i)+h*l1)-y*(I(i)+h*l1);
m2=y*(I(i)+h*l1);