MATLAB: I need help solving an ODE using the Runge Kutta method. The values I obtain for “A” are correct, but the value I receive for “C” are incorrect.

homeworkoderunge kutta

% Constants
dt = 1; % Step size in hours
tf = 50; % Solve from time 0 to 50
Ka = 0.20; % rate of absorption at 1/hr
Ke = 0.047; % rate of elimination at 1/hr
Vd = 50; % volume of distribution in mL
% Define the ODE
f1 = @(A,t) -Ka*A;
f2 = @(A,C,t) Ka*A/Vd-(Ke*C);
% Initial Condition
A(1) = 250; % Initial Dose of Fenitoine in µg
t(1) = 0;
C(1) = 0;
% RK for loop
 for k = 1:tf;
    t(k+1) = t(k)+dt; 
    k1(k) = dt*f1(A(k),t(k)); 
    k2(k) = dt*f1((A(k)+0.5*k1(k)),t(k)+0.5*dt);    
    k3(k) = dt*f1((A(k)+0.5*k2(k)),t(k)+0.5*dt); 
    k4(k) = dt*f1((A(k)+k3(k)),t(k)+dt); %    
    A(k+1) = A(k)+(k1(k)+ 2*k2(k)+ 2*k3(k)+ k4(k))/6; 
 end
 A; % Mass of drug in ?g remaining to be absorbed in a period of 50 hours
for k = 1:tf;
    t(k+1) = t(k)+dt;
    k5(k) = dt*f2(A(k),C(k),t(k)); 
    k6(k) = dt*f2((A(k)+0.5*k1(k)),t(k)+0.5*dt,C(k)+0.5*k5(k));
    k7(k) = dt*f2((A(k)+0.5*k2(k)),t(k)+0.5*dt,C(k)+0.5*k6(k)); 
    k8(k) = dt*f2((A(k)+k3(k)),t(k)+dt,C(k)+k7(k)); 
    C(k+1) = C(k)+(k5(k)+ 2*k6(k)+ 2*k7(k)+ k8(k))/6; 
end
C; % Drug concentration
These are the values I should be getting for C
 0 
  0.8847 
  1.5684 
  2.0894 
  2.4790 
  2.7627 
  2.9613 
  3.0918 
  3.1680 
  3.2011 
  3.2004 
  3.1732 
  3.1255 
  3.0623 
  2.9874 
  2.9040 
  2.8148 
  2.7216 
  2.6262 
  2.5298 
  2.4334 
  2.3379 
  2.2438 
  2.1517 
  2.0618 
  1.9744 
  1.8897 
  1.8078 
  1.7288 
  1.6527 
  1.5795 
  1.5092 
  1.4417 
  1.3770 
  1.3149 
  1.2556 
  1.1987 
  1.1443 
  1.0923 
  1.0426 
  0.9951 
  0.9497 
  0.9064 
  0.8650 
  0.8254 
  0.7876 
  0.7516 
  0.7172 
  0.6843 
  0.6530 
  0.6230

Best Answer

Your signature of the derivative function is f2(A,C,t) not f2(A,t,C). So these lines:
    k6(k) = dt*f2((A(k)+0.5*k1(k)),t(k)+0.5*dt,C(k)+0.5*k5(k));
    k7(k) = dt*f2((A(k)+0.5*k2(k)),t(k)+0.5*dt,C(k)+0.5*k6(k)); 
    k8(k) = dt*f2((A(k)+k3(k)),t(k)+dt,C(k)+k7(k));
should probably be these instead:
    k6(k) = dt*f2((A(k)+0.5*k1(k)),C(k)+0.5*k5(k),t(k)+0.5*dt);
    k7(k) = dt*f2((A(k)+0.5*k2(k)),C(k)+0.5*k6(k),t(k)+0.5*dt); 
    k8(k) = dt*f2((A(k)+    k3(k)),C(k)+    k7(k),t(k)+    dt);
Related Question