ОДУ СИСТЕМА УРАВНЕНИЯ

Question

Я пытаюсь решить систему од, используя Runge-kutta, я сделал функцию для RK2(f,h,x0,y0,xfinal) и попытался решить систему, показанную ниже, с указанными I C. Может ли кто-нибудь помочь исправить код, поскольку я получаю ошибки и код не работает. оды

beta = 1/3;
gamma = 1/7;
syms R S I % Symbolic Math Toolbox
N = S+I+R;
ode1 = -(beta*I*S)/N;
ode2 = -(beta*I*S)/N-gamma*I;
ode3 = gamma*I;
odes = [ode1,ode2,ode3];
for j = odes
    RK2(j,0.2,0,8e6,7);
end


function [xs,ys] = RK2(f,h,x0,y0,xfinal)
ffnc = matlabFunction(f);
fprintf('\n x        y ');
o = 1;
while x0 <= xfinal
    fprintf('\n%4.3f  %4.3f ',x0,y0); %values of x and y
    xs(o) = x0;
    ys(o) = y0;
    k1 = h*ffnc (x0,y0); 
    x1 = x0+h; 
    k2 = h*ffnc (x1,y0+k1);
    y1 = y0+(k1+k2)/2;           
    x0 = x1;
    y0 = y1;  
    o = o+1;
end
end