Modeling the radio-active decay using ode23

Hi I got a question like this. The basic equation for modeling radio-active decay is : dx/dt =−rx where x is the amount of the radio-active substance at time t and r is the decay rate. Some radio-active substances decay into other radioactive substances which in turn also decay. For example, Strontium 92 (r1 = 0.256/h) decays into Yttrium 92 (r2 = 0.127/h), which in turn decays in to Zirconium. Write down a pair of differential equations for Strontium and Yttrium to describe what is happening. Starting at t=0 with 5 x1026 atoms of Strontium 92 and none of Yttrium, use the Runge-Kutta method (ode23) to solve the equations up to t=8 in steps of 1/3 h. Plot the results.

So I'm trying to solve it simultaneously but obviously the code has some bug. I do not know how to subtract decayed weigh from the original rate. The error probably happens in dxdt function. Could anyone help me out? Thank you.

 function matlab
 clc;clear;

 %Radioactive decay
 y0=[5*10^26;0];
 soln = ode23(@f1,[0 8],y0)
 t = linspace(0,8,24);
 y(:,1)=deval(soln,t,1); %Strontium
 y(:,2)=deval(soln,t,2); % Yttrium

 figure
 plot(t,y(:,1),'-o',t,y(:,2),'--');
 hold on;grid on;
 legend('Strontium','Yttrium')
 end

 function dxdt = f1(x,t)
 r1 = 0.256;
 r2 = 0.127;
 dxdt(1) = -r1 * x;
 dxdt(2) = -r2 * x;
 dxdt =dxdt';
 end

Answer 1

[James Tursa]

For starters, you have t and x reversed in the argument list of the derivative function f1. Also, it seems to me that the rate of change of Yttrium should be a combination of how much current Yttrium is decaying (a Yttrium loss) and how much Strontium is decaying into Yttrium (a Yttrium gain). So tweaking your code a bit, maybe something like this is what you need:

   function raddecay
   clc;clear;

   %Radioactive decay
   y0=[5*10^26;0];
   t = [0 8];
   [T,Y] = ode23(@f1,t,y0);

   figure
   plot(T,Y(:,1),'-o',T,Y(:,2),'--');
   hold on;grid on;
   legend('Strontium','Yttrium')
   end

   function dxdt = f1(t,x)
   dxdt = [0;0];
   r1 = 0.256;
   r2 = 0.127;
   dxdt(1) = -r1 * x(1);  % Strontium decay
   dxdt(2) = -r2 * x(2) + r1 * x(1); % Yttrium decay + gain from Strontium decay
   end