How to obtain the numerical solution of these differential equations with matlab

I have differential equations derived from epidemic spreading. I want to obtain the numerical solutions. Here's the equations,

t is a independent variable and ranges from [0,100]. The initial value is

y1 = 0.99; y2 = 0.01; y3 = 0;

At first, I planned to deal these with ode45 function in matlab, however, I don't know how to express the series and the combination. So I'm asking for help here.

**

The problem is how to express the rig开发者_运维知识库ht side of the equations as the odefun, which is a parameter in the ode45 function.

**

Matlab has functions to calculate binomial coefficients (number of combinations) and the finite series can be expressed just as matrix multiplication. I'll demonstrate how that works for the sum in the first equation. Note the use of the element-wise "dotted" forms of the arithmetic operators.

Calculate a row vector coefs with the constant coefficients in the sum as:

octave-3.0.0:33> a = 0:20;
octave-3.0.0:34> coefs = log2(a * 0.05 + 1) .* bincoeff(20, a);

The variables get combined into another vector:

octave-3.0.0:35> y1 = 0.99;
octave-3.0.0:36> y2 = 0.01;
octave-3.0.0:37> z = (y2 .^ a) .* ((1 - y2) .^ a) .* (y1 .^ a);

And the sum is then just evaluated as the inner product:

octave-3.0.0:38> coefs * z'

The other sums are similar.

 function demo(a_in)
 X = [0;0;0];
  T = [0:.1:100];
  a = a_in; % for nested scope

  [Xout,  Tout ]= ode45( @myFunc, T, X );

    function [dxdt] = myFunc( t, x )
         % nested function accesses "a"
         dxdt = 0*x + a; 
         % Todo: real value of dxdt. 
    end
  end

What about this, and you simply need to fill in the dxdt from your math above? It remains to be seen if the numerical roundoff matters...

Edit: there's a serious issue due to the 1=y1+y2+y3 constraint. Is that even allowed, since you have an IVP with 3 initial values given and 3 first order ODE's? If that constraint is a natural consequence of the equations, it may not be needed.