ODE45 Solution for a coupled equation

Hi, I have a coupled differential eqaution for a decay model. I have tried to solve the problem with ODE45 (code given below). The solution of this eqaution S and T should follow (S+T) ≤ (S0+T0) at any point of tspan. Though I am not getting a solution like that. Please let me know why and how can fix this. Thanks.

function ODEDecayModel()
clear all, close all,  clc
  function dotn = Mat (tspan,n)
      k_rs = 1.07e-02; 
     k_ISC = 1.33e-02; 
     k_RISC = 1.13e-03; 
     k_nT = 3.53e-4; 
 dotn = zeros(2,1);
     dotn(1) = -k_rs*n(1)-k_ISC*n(1)+k_RISC*n(2);
     dotn(2) = -k_RISC*n(2)+k_ISC*n(1)-k_nT*n(2);
  end

tspan = 00:0.1:1800;
S0 = 5.80e18; % S at 0
T0 = 0.1;% T at 0
IC = [S0 T0];%intial condition

[t,Values] = ode45(@Mat,tspan,IC); 
S = Values(:,1)
T = Values(:,2)

figure(1)
plot(t,S,'b')
hold on
plot(t,T,'r')
set(gca,'xscale','log')
end

[Bjorn Gustavsson]

What purpose does your clear all, close all clc serve? Inside a function taking no input arguments...

[Monirul Hasan]

Didn't noticed it at all. Thanks

Answer 1

[Bjorn Gustavsson]

When running your function I get an initial growth of T, which is due to the k_ISC*n(1) term in your ode - that would represent decay of S to T, then after an initial growth of T (compare the initial condition for S and T!) it starts to decay too (with some contribution back to S). If you change your plot to:

plot(t,[Values,sum(Values')'])
set(gca,'yscale','log')
axis([0 1800 1e16 1e19])

You will see that the sum of S and T does indeed decrease at all times.

HTH

[Monirul Hasan]

Dear Bjorn,

Thanks a lot. I get it.