runge kutta method 4th order

clc;                                               % Clears the screen
clear all;
lamda=0.2;
mu=0.8;
h=0.1;                                             % step size
x = 0:h:1;                                         % Calculates upto y(1)
y = zeros(1,length(x)); 
z = zeros(1,length(x));
% It calculates ODE using Runge-Kutta 4th order method
% Author Ido Schwartz
% Originally available form: http://www.mathworks.com/matlabcentral/fileexchange/29851-runge-kutta-4th-order-ode/content/Runge_Kutta_4.m
% Edited by Amin A. Mohammed, for 2 ODEs(April 2016)
clc;                                               % Clears the screen
clear all;
h=0.1;                                             % step size
x = 0:h:1;                                         % Calculates upto y(1)
y = zeros(1,length(x)); 
z = zeros(1,length(x)); 
y(1) = 3;                                          % initial condition
z(1) = 1; 
% It calculates ODE using Runge-Kutta 4th order method
% Author Ido Schwartz
% Originally available form: http://www.mathworks.com/matlabcentral/fileexchange/29851-runge-kutta-4th-order-ode/content/Runge_Kutta_4.m
% Edited by Amin A. Mohammed, for 2 ODEs(April 2016)
clc;                                               % Clears the screen
clear all;
h=0.1;                                             % step size
x = 0:h:1;                                         % Calculates upto y(1)
y = zeros(1,length(x)); 
z = zeros(1,length(x)); 
y(1) = 0.2;                                        % initial condition
z(1) = 0.3; 
m(1)=  0.2;
% initial condition
%F_xyz = @(x,y,z)  (lamda^2/(lamda+mu)^2)*exp(-2*(lamda+mu)*t)+((2*mu*lamda)/(lamda+mu)^2)*exp(-(mu+lamda)*t)+mu^2/(lamda+mu)^2;                                  % change the function as you desire
%G_xyz = @(x,y,z) ((2*mu*lamda)/(lamda+mu)^2)+(((2*lamda*(lamda-mu))/(lamda+mu)^2)*exp(-(mu+lamda)*t))-2*((lamda^2/lamda+mu)^2)*exp(-2*(mu+lamda)*t);
%O_xyx=@(x',y,z) (lamda^2/(lamda+mu)^2)*exp(-2*(mu+lamda)*t))-(2*lamda^2/(lamda+mu)^2)*exp(-(mu+lamda)*t)+lamda^2/(lamda+mu)^2;
for i=1:(length(x)-1)                              % calculation loop
    k_1 = F_xyz(x(i),y(i),z(i));
    L_1 = G_xyz(x(i),y(i),z(i));
    k_2 = F_xyz(x(i)+0.5*h,y(i)+0.5*h*k_1,z(i)+0.5*h*L_1);
    L_2 = G_xyz(x(i)+0.5*h,y(i)+0.5*h*k_1,z(i)+0.5*h*L_1);
    k_3 = F_xyz((x(i)+0.5*h),(y(i)+0.5*h*k_2),(z(i)+0.5*h*L_2));
    L_3 = G_xyz((x(i)+0.5*h),(y(i)+0.5*h*k_2),(z(i)+0.5*h*L_2));
    k_4 = F_xyz((x(i)+h),(y(i)+k_3*h),(z(i)+L_3*h)); % Corrected        
    L_4 = G_xyz((x(i)+h),(y(i)+k_3*h),(z(i)+L_3*h));
    y(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;  % main equation
    z(i+1) = z(i) + (1/6)*(L_1+2*L_2+2*L_3+L_4)*h;  % main equation
end

if i have 3 initial value (0.2,0.3,0.2) and for each value there is an equation

F_xyz for first; G_xyz for second and O_xyz for third

how can i merge the 3 value in one plot?

i get confused to solve it

Answer 1

[Sulaymon Eshkabilov]

Hi,

There are several flaws and unclear issues in your script:

(1) What is your given system of ODEs with x,y,z variables?

(2) F_xyz, G_xyz, O_xyz need to be defined properly via anonymous function, function file, inline function or directly or matlabFunction. They are not defined properly within your script.

(3) within the loop [for .. end], x, y, z need to be calculated. In your script, x is missing and O_xyz is missing. From your statements, presumably you are trying to compute x, y, z. The defined x =0:h:1 is time space???

(4) ICs for x0, y0, z0 need to be defined correctly

(5) Memory allocation can be adjusted a bit more accurately, e.g. y = [y(1), zeros(1, numel(x)-1)]; ...

Good luck.