ODE45 solving 3 equations and comparing to RK-4

%% RK4 methods
%% Description
% RK4 methods are used to solve coupled first-order ODEs
clear
close all
clc

%% Input variables
y10 = 10; % initial value of T*
y20 = 100; % initial value of V1
y30 = 0; % initial value of Vx
tf = 5; % interval size (day)
h = 1.0; % step size (day)
t = 0:h:tf; % vector of time points
numPoints = length(t);

IC_vector = [ y10 ; y20 ; y30] ; 

%% Calculations
%apply initial conditions
y1(1) = y10; % initial condition at t = 0
y2(1) = y20; % initial condition at t = 0
y3(1) = y30; % initial condition at t = 0


% RK-4 method for solving coupled first-order ODEs 
for k = 1:numPoints-1

    [k11, k21, k31] = func_hiv(t(k), y1(k), y2(k), y3(k)); % <--- your input
    [k12, k22, k32] = func_hiv( t(k)+ h/2, y1(k)+0.5*h*k11 , y2(k)+0.5*h*k21 , y3(k)+0.5*h*k31 ); 
    [k13, k23, k33] = func_hiv( t(k)+ h/2 , y1(k)+0.5*h*k12 , y2(k)+0.5*h*k22 , y3(k)+ 0.5*h*k32 );
    [k14, k24, k34] = func_hiv( t(k) , y1(k)+0.5*h*k13 , y2(k)+ 0.5*h*k23, y3(k)+0.5*h*k33 );

    k1avg = (1/6)*(k11+k12+k13+k14 ); % <--- your input
    k2avg = (1/6) * (k21+k22+k23+24 );
    k3avg = (1/6)* (k31+k32+k33+k34 );

    y1(k+1) = y1(k) + k1avg*(t(k+1)-t(k)); % <--- your input
    y2(k+1) = y2(k) + k2avg*(t(k+1)-t(k));
    y3(k+1) = y3(k) + k3avg*(t(k+1)-t(k));

end

% ODE45 built-in function
% <--- your input

% ode_sys = [( k*y2*Tc-y1*delta ) ; (-c*y2) ; (N*delta*y1-c*y3   ) ] ; 
 [t,y] = ode45(@(t,y) func_hiv(t,y1,y2,y3),t,IC_vector) ; 

Error using odearguments
@(T,Y)FUNC_HIV(T,Y1,Y2,Y3) must return a column vector.

Error in ode45 (line 104)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

% Plot Results
subplot(1,3,1)
plot(t,y1,'k-','LineWidth',2)
hold on 
 % <--- your input
set(gca,'FontSize',16,'LineWidth',2)
xlabel('{\itt} in day')
ylabel('{\itT}^* infected cells/ \muL')
legend('RK-4', 'ode45');

subplot(1,3,2)
plot(t,y2,'k-','LineWidth',2)
hold on 
 % <--- your input 
set(gca,'FontSize',16,'LineWidth',2)
xlabel('{\itt} in day')
ylabel('{\itV}_1 virions/ \muL')
legend('RK-4', 'ode45');


subplot(1,3,3)
plot(t,y3,'k-','LineWidth',2)
hold on
 % <--- your input
set(gca,'FontSize',16,'LineWidth',2)
xlabel('{\itt} in day')
ylabel('{\itV}_X noninfectious virus particles/ \muL')
legend('RK-4', 'ode45');


%% Functions
function [f1, f2, f3] = func_hiv(t,  y1, y2, y3)
    % Evaluate slopes of the coupled ODEs for HIV-1 dynamics in bloodstream
    % Constants
    k = 2e-4; % uL/day/virions
    N = 60; % virions/cell
    delta = 0.5; % /day
    c = 3.0; % /day
    Tc = 250; % noninfected cells/uL

    % Equations % <--- your input
     
    f1 =  k*y2*Tc-y1*delta ; 
    f2 = (-c*y2); 
    f3 =   (N*delta*y1-c*y3   ) ;

    
end

I have tried creating a seperate function that uses dydt as an output that containts each of the equations , and it gives me an error due to it counting each constant in the vector, not maching the 3 initial conditions, we are supposed to use the given function which defines [f1 f2 f3] . The code runs fine without the ode45 but it is needed to compare, thanks !

Answer 1

[Torsten]

%% RK4 methods
%% Description
% RK4 methods are used to solve coupled first-order ODEs
clear
close all
clc

%% Input variables
y10 = 10; % initial value of T*
y20 = 100; % initial value of V1
y30 = 0; % initial value of Vx
tf = 5; % interval size (day)
h = 1.0; % step size (day)
t = 0:h:tf; % vector of time points
numPoints = length(t);

IC_vector = [ y10 ; y20 ; y30] ; 

%% Calculations
%apply initial conditions
y1(1) = y10; % initial condition at t = 0
y2(1) = y20; % initial condition at t = 0
y3(1) = y30; % initial condition at t = 0


% RK-4 method for solving coupled first-order ODEs 
for k = 1:numPoints-1

    [k11, k21, k31] = func_hiv(t(k), y1(k), y2(k), y3(k)); % <--- your input
    [k12, k22, k32] = func_hiv( t(k)+ h/2, y1(k)+0.5*h*k11 , y2(k)+0.5*h*k21 , y3(k)+0.5*h*k31 ); 
    [k13, k23, k33] = func_hiv( t(k)+ h/2 , y1(k)+0.5*h*k12 , y2(k)+0.5*h*k22 , y3(k)+ 0.5*h*k32 );
    [k14, k24, k34] = func_hiv( t(k) , y1(k)+0.5*h*k13 , y2(k)+ 0.5*h*k23, y3(k)+0.5*h*k33 );

    k1avg = (1/6)*(k11+k12+k13+k14 ); % <--- your input
    k2avg = (1/6) * (k21+k22+k23+24 );
    k3avg = (1/6)* (k31+k32+k33+k34 );

    y1(k+1) = y1(k) + k1avg*(t(k+1)-t(k)); % <--- your input
    y2(k+1) = y2(k) + k2avg*(t(k+1)-t(k));
    y3(k+1) = y3(k) + k3avg*(t(k+1)-t(k));

end

% ODE45 built-in function
% <--- your input

% ode_sys = [( k*y2*Tc-y1*delta ) ; (-c*y2) ; (N*delta*y1-c*y3   ) ] ; 
 [t,y] = ode45(@func_hiv_help,t,IC_vector) ; 


% Plot Results
subplot(1,3,1)
plot(t,y1,'k-','LineWidth',2)
hold on 
 % <--- your input
set(gca,'FontSize',16,'LineWidth',2)
xlabel('{\itt} in day')
ylabel('{\itT}^* infected cells/ \muL')
legend('RK-4', 'ode45');

Warning: Ignoring extra legend entries.

subplot(1,3,2)
plot(t,y2,'k-','LineWidth',2)
hold on 
 % <--- your input 
set(gca,'FontSize',16,'LineWidth',2)
xlabel('{\itt} in day')
ylabel('{\itV}_1 virions/ \muL')
legend('RK-4', 'ode45');

Warning: Ignoring extra legend entries.

subplot(1,3,3)
plot(t,y3,'k-','LineWidth',2)
hold on
 % <--- your input
set(gca,'FontSize',16,'LineWidth',2)
xlabel('{\itt} in day')
ylabel('{\itV}_X noninfectious virus particles/ \muL')
legend('RK-4', 'ode45');

Warning: Ignoring extra legend entries.

%% Functions
function [f1, f2, f3] = func_hiv(t,  y1, y2, y3)
    % Evaluate slopes of the coupled ODEs for HIV-1 dynamics in bloodstream
    % Constants
    k = 2e-4; % uL/day/virions
    N = 60; % virions/cell
    delta = 0.5; % /day
    c = 3.0; % /day
    Tc = 250; % noninfected cells/uL

    % Equations % <--- your input
     
    f1 =  k*y2*Tc-y1*delta ; 
    f2 = (-c*y2); 
    f3 =   (N*delta*y1-c*y3   ) ;

    
end

function f = func_hiv_help(t,y)
  [f1,f2,f3] = func_hiv(t,  y(1), y(2), y(3));
  f = [f1;f2;f3];
end

And k1avg, k2avg and k3avg are set wrong. Look up again the RK-4 scheme to see how they are computed.