I dont know what to do

% Fourth Order Runge Kutta Method to solve system of ODEs
%  y1'= f1(t, y1, y2, y3, y4);
%  y2'= f2(t, y1, y2, y3, y4);
%  y3'= f2(t, y1, y2, y3, y4);
%  y4'= f2(t, y1, y2, y3, y4);

% Fourth Order RK method for solving ODE
% y' = f(t, y)
% y_n+1 =y_n + (h/6) * (delta_y1 + 2*delta_y2 + 2*delta_y3 + delta_y4)
% delta_y1 = f(tn,yn)
% delta_y2 = f(tn + h/2, yn + h*(delta_y1/2))
% delta_y3 = f(tn + h/2, yn + h*(delta_y2/2))
% delta_y4 = f(tn + h, yn + h*delta_y3)

% dy1/dx = -*y1*(y2+y3)
% dy2/dx = 0.47*(1.41*10^-9)*y1*(y2+y3)-0.1*y2 
% dy3/dx = 0.53*(1.41*10^-9)*y1*(y2+y3)-0.1*y3
% dy4/dx = 4 - 0.3 * y2 - 0.1 * y1

% Numerically integrate above from x = 0 to x = 60 using step size of 0.001
% IC: y1(x=0) = 8.18*10**7; y2(x=0) = 5.0604*10**3; y3(x=0) = 4.4876*10**3; y4(x=0) = 4590;
% y1exact = ;
% y2exact = ;
% y3exact = ;
% y4exact = ;

clear
clc
clf

% Initial Value and time step
a = 1.41*10^-9; % Lower limit of Integration
b = 2*10^-9; % Upper limit of Integration
h = 0.001; % h = delta_x
y1(1) = 8.18*10^7; % S(1)
y2(1) = 5.0604*10^3; % I_1(1)
y3(1) = 4.4876*10^3; % I_2(2)
y4(1) = 4590; % R(1)
n = (b-a)/h + 1; % no. of points

% Solution
x(1) = 0;
func1 = @(y1, y2, y3, y4) (-1.41*10^-9 * y1 * (y2 + y3)); % function to integrate dS(t)/dx
func2 = @(y1, y2, y3, y4) ((0.45 * 1.41*10^-9 * y1 * (y2 + y3)) - ((1/10) * y2)); % function to integrate dI_1(t)/dx
func3 = @(y1, y2, y3, y4) (((1 - 0.45) * 1.41*10^-9 * y1 * (y2 + y3)) - (((1/15) + (1/50)) * y3)); % function to integrate dI_2(t)/dx
func4 = @(y1, y2, y3, y4) (((1/10) * y2) + ((1/15) * y3)); % function to integrate dR(t)/dt
% func_exact = @(x) ;

for i = 1:n-1
    x(i+1) = 1*h;
    
    delta_y1_1 = h * feval(func1, y1(i));
    delta_y2_1 = h * feval(func1, y1(i) + (delta_y1_1)/2);
    delta_y3_1 = h * feval(func1, y1(i) + (delta_y2_1)/2);
    delta_y4_1 = h + feval(func1, y1(i) + delta_y3_1);
    y1(i+1) = y1(i) + (1/6) + (delta_y1_1 + 2*delta_y2_1 + 2*delta_y3_1 + delta_y4_1);
    
    delta_y1_2 = h * feval(func2, y1(i), y2(i));
    delta_y2_2 = h * feval(func2, y1(i) + (delta_y1_1)/2, y2(i) + (delta_y1_2)/2);
    delta_y3_2 = h * feval(func2, y1(i) + (delta_y2_1)/2, y2(i) + (delta_y2_2)/2);
    delta_y4_2 = h + feval(func2, y1(i) + delta_y3_1, y2(i) + delta_y3_2);
    y2(i+1) = y2(i) + (1/6) + (delta_y1_2 + 2*delta_y2_2 + 2*delta_y3_2 + delta_y4_2);
    
    delta_y1_3 = h * feval(func3, y1(i), y2(i), y3(i));
    delta_y2_3 = h * feval(func3, y1(i) + (delta_y1_1)/2, y2(i) + (delta_y1_2)/2, y3(i) + (delta_y1_3)/2);
    delta_y3_3 = h * feval(func3, y1(i) + (delta_y2_1)/2, y2(i) + (delta_y2_2)/2, y3(i) + (delta_y2_3)/2);
    delta_y4_3 = h + feval(func3, y1(i) + delta_y3_1, y2(i) + delta_y3_2, y3(i) + delta_y3_3);
    y3(i+1) = y3(i) + (1/6) + (delta_y1_3 + 2*delta_y2_3 + 2*delta_y3_3 + delta_y4_3);

    delta_y1_4 = h * feval(func4, y1(i), y2(i), y3(i), y4(i));
    delta_y2_4 = h * feval(func4, y1(i) + (delta_y1_1)/2, y2(i) + (delta_y1_2)/2, y3(i) + (delta_y1_3)/2, y4(i) + (delta_y1_4)/2);
    delta_y3_4 = h * feval(func4, y1(i) + (delta_y2_1)/2, y2(i) + (delta_y2_2)/2, y3(i) + (delta_y2_3)/2, y4(i) + (delta_y2_4)/2);
    delta_y4_4 = h + feval(func4, y1(i) + delta_y3_1, y2(i) + delta_y3_2, y3(i) + delta_y3_3, y4(i) + delta_y3_4);
    y4(i+1) = y4(i) + (1/6) + (delta_y1_4 + 2*delta_y2_4 + 2*delta_y3_4 + delta_y4_4);

end
x;
y1;
y2;
y3;
y4;

[Walter Roberson]

What problems are you having?

[SALAH ALRABEEI]

what to do

[Meor Hasan Meor Jumat]

the results doesnt come out

Answer 1

[Jan]

Remove the semicolons from these lines:

x;
y1;
y2;
y3;
y4;

The trailing semicolon suppresses the output. Therefore "x;" shows x without showing x. Of course you do not see anything. You will see, that all variables are scalars, because you calculate 0 steps

a = 1.41*10^-9; % Lower limit of Integration
b = 2*10^-9;    % Upper limit of Integration
h = 0.001;      % h = delta_x
n = (b - a) / h + 1

Now n = 1.00000059

Then:

for i = 1:n-1

performs no iterations and the trajectory stays at the initial values.

By the way: -1.41*10^9 is an expensive power operations, while -1.41e9 is a cheap constant.

I'd prefer to convert:

    func4 = @(y1, y2, y3, y4) (((1/10) * y2) + ((1/15) * y3));
    % to
    func4 = @(y1, y2, y3, y4) y2 / 10 + y3 / 15;    

There is something wrong with your Runge-Kutta code. Why does the 2nd component depend on the result of the 1st one, but not the other way around?

func1 = @(y1, y2, y3, y4), so it depends on 4 inputs. But here you provide only 1 input:

    delta_y1_1 = h * feval(func1, y1(i));    

This bug did not cause problems yet, because the loop is not entered at all.

Calculating the functions elementwise demands for writing a Runge Kutta code for each number of dimensions. It is much easier to use vectors instead:

y = zeros(n, 4);
x = zeros(n, 1);
y(1, :) = [8.18e7, 5.0604e3; 4.4876e3; 4590];
x(1)    = a;  % Not 0

fcn = @(x, y) [-1.41e-9 * y(1) * (y(2) + y(3)), ...
               0.45 * 1.41e-9 * y(1) * (y(2) + y(3)) - y(2) / 10, ...
               0.55 * 1.41e-9 * y(1) * (y(2) + y(3)) - y(3) * 13 / 150, ...
               y(2) / 10 + y(3) / 15];

for i = 1:n-1
   x(i + 1) = x(i) + h;   % Not: "x(i+1) = 1*h", but maybe: "i * h"
   d1 = h * fcn(x(i), y(i, :));
   d2 = h * fcn(x(i), y(i, :) + d1 / 2);
   d3 = h * fcn(x(i), y(i, :) + d2 / 2);
   d4 = h + fcn(x(i), y(i, :) + d2);
   y1(i+1, :) = y1(i, :) + (d1 + d2 + d2 + d4) / 6;
end

Do you see the benefit of using vectors? Less clutter, less bugs.

Although your fcn does not depened on x, I've inserted it here such that you can use the Runge Kutta part in general.

[Meor Hasan Meor Jumat]

Thank you for respond. I did the correction but still have errors. It's said "Error using zeros

Size inputs must be integers.

Error in line 49

y = zeros(n, 4)"

Can you give me a feedback.

[Image Analyst]

So your n = 1.00000059. What were you hoping it would be? 1? Something else?

[Walter Roberson]

n = (b-a)/h + 1; % no. of points

Your current calculation gives you a floating point non-integer number for n. But MATLAB cannot have (for example) 1.000009 array locations. You need to decide whether you want to round the calculated value down with floor() or up with ceil() .

[James Tursa]

@Meor Hasan Meor Jumat Please post your current code. All variables need to be updated simultaneously, so I would suggest you take Jan's advice and use state vectors to accomplish this.

[Meor Hasan Meor Jumat]

clear
clc
clf

% Initial Value and time step
a = 21; % Lower limit of Integration
b = 31; % Upper limit of Integration
h = 0.5; % h = delta_x
n = (b-a)/h + 1; % no. of points

y = zeros(n, 4);
x = zeros(n, 1);
y(1, :) = [10276; 751; 11; 9395];
x(1) = a; % Not 0

fcn = @(x, y) [-10219 * y(1) * (y(2) + y(3)),...
                1.26 * 10219 * y(1) * (y(2) + y(3)) - 7.31 * y(2), ...
                -0.26 * 10219 * y(1) * (y(2) + y(3)) - (0.11 + 1.26) * y(3), ...
                7.31 * y(2) + 0.11 * y(3) ]; % system function

for i = 1:n-1
    x( i + 1 ) = x(i) + h; % Not: "x(i+1) = 1*h", but maybe: "i * h"
    
    d1 = h * fcn(x(i), y(i, :));
    d2 = h * fcn(x(i), y(i, :) + d1 / 2);
    d3 = h * fcn(x(i), y(i, :) + d2 / 2);
    d4 = h + fcn(x(i), y(i, :) + d2);
    y(i+1, :) = y(i, :)  + (d1 + d2 + d3 + d4) / 6;
    
end
x
y
plot ( x, y )

[Walter Roberson]

For demonstration, I set an upper limit of n = 5.

Notice the exponents of the numbers in y by the time you get to 5 steps.

% Initial Value and time step
a = 21; % Lower limit of Integration
b = 31; % Upper limit of Integration
h = 0.5; % h = delta_x
n = (b-a)/h + 1; % no. of points
n = 5;

y = zeros(n, 4, 'sym');
x = zeros(n, 1, 'sym');
y(1, :) = [10276; 751; 11; 9395];
x(1) = a; % Not 0

fcn = @(x, y) [-10219 * y(1) * (y(2) + y(3)),...
                1.26 * 10219 * y(1) * (y(2) + y(3)) - 7.31 * y(2), ...
                -0.26 * 10219 * y(1) * (y(2) + y(3)) - (0.11 + 1.26) * y(3), ...
                7.31 * y(2) + 0.11 * y(3) ]; % system function

for i = 1:n-1
    x( i + 1 ) = x(i) + h; % Not: "x(i+1) = 1*h", but maybe: "i * h"
    
    d1 = h * fcn(x(i), y(i, :));
    d2 = h * fcn(x(i), y(i, :) + d1 / 2);
    d3 = h * fcn(x(i), y(i, :) + d2 / 2);
    d4 = h + fcn(x(i), y(i, :) + d2);
    y(i+1, :) = y(i, :)  + (d1 + d2 + d3 + d4) / 6;
    
end
x

x = 

y

y = 

vpa(y)

ans = 

plot ( x, y )

[James Tursa]

This

y(i+1, :) = y(i, :)  + (d1 + d2 + d3 + d4) / 6;

should be this

y(i+1, :) = y(i, :)  + (d1 + 2*d2 + 2*d3 + d4) / 6;

[Walter Roberson]

Using the hint from James, let us experiment:

% Initial Value and time step
a = 21; % Lower limit of Integration
b = 31; % Upper limit of Integration
h = 0.5; % h = delta_x
n = (b-a)/h + 1; % no. of points
n = 5;

y = zeros(n, 4, 'sym');
x = zeros(n, 1, 'sym');
y(1, :) = [10276; 751; 11; 9395];
x(1) = a; % Not 0

fcn = @(x, y) [-10219 * y(1) * (y(2) + y(3)),...
                1.26 * 10219 * y(1) * (y(2) + y(3)) - 7.31 * y(2), ...
                -0.26 * 10219 * y(1) * (y(2) + y(3)) - (0.11 + 1.26) * y(3), ...
                7.31 * y(2) + 0.11 * y(3) ]; % system function

for i = 1:n-1
    x( i + 1 ) = x(i) + h; % Not: "x(i+1) = 1*h", but maybe: "i * h"
    
    d1 = h * fcn(x(i), y(i, :));
    d2 = h * fcn(x(i), y(i, :) + d1 / 2);
    d3 = h * fcn(x(i), y(i, :) + d2 / 2);
    d4 = h + fcn(x(i), y(i, :) + d2);
    y(i+1, :) = y(i, :)  + (d1 + 2*d2 + 2*d3 + d4) / 6;
    
end
x

x = 

y

y = 

vpa(y)

ans = 

plot ( x, y )

This suggests that there are other problems with the system.