Solving ODE with big dimension

I have an ODE code to run an SEIR type epidemiological system segemented by age (n=21), and social group (soc=10). I expect the output to be of the form Classes.S(:,:,i) where i=1:10. The code is runing but only the first row of my initail condition is repeated 10times.

This is the code snippet:

function [Classes] = ODE_social_TEST2(para,ICs,maxtime,Mix_H,Mix_W,Mix_S,Mix_O,n,soc,w)

     opts = odeset('RelTol',1e-4,'AbsTol',1e-3); %tolerance level%tolerance level
     tspan = [0:1:maxtime]; %tie span
      
     %Shape of the initial condition

     All = reshape([ICs.S'; ICs.E'; ICs.A'; ICs.J';ICs.Q';ICs.I';ICs.H'; ICs.R'; ICs.D'; ICs.Ir'], []*n,1);
    
    
     [t , pop] = ode15s(@(t,y)diff_socialmodel(t,y,para),tspan,All,opts);

            
     %% Define the structure of the output

     Classes = struct('S',zeros(numel(t), n, 1),'E',zeros(numel(t), n, 1),'A',zeros(numel(t), n, 1),'J',zeros(numel(t), n, 1),'Q',zeros(numel(t), n, 1), ...
         'I',zeros(numel(t), n, 1),'H',zeros(numel(t), n, 1),'R',zeros(numel(t), n, 1),'D',zeros(numel(t), n, 1),'Ir',zeros(numel(t), n, 1),'t',[]);

 

     for i = 1:soc
         Classes.S(:,:,i) = pop(:,1:1:n);
         Classes.E(:,:,i) = pop(:,n+1:1:2*n);
         Classes.A(:,:,i) = pop(:,(2*n)+1:1:3*n);
         Classes.J(:,:,i) = pop(:,(3*n)+1:1:4*n);
         Classes.Q(:,:,i) = pop(:,(4*n)+1:1:5*n);
         Classes.I(:,:,i) = pop(:,(5*n)+1:1:6*n);
         Classes.H(:,:,i) = pop(:,(6*n)+1:1:7*n);
         Classes.R(:,:,i) = pop(:,(7*n)+1:1:8*n);
         Classes.D(:,:,i) = pop(:,(8*n)+1:1:9*n);
         Classes.Ir(:,:,i) = pop(:,(9*n)+1:1:10*n);
     end

     Classes.t = t;

     %% Set up the population change

    function dpop = diff_socialmodel(t,pop,para)

       
        %% DEfine the age structured mixing

        ageMix = Mix_H + Mix_W + Mix_S + Mix_O;


       
        S = pop(1:1:n,:);
        E = pop(n+1:1:2*n,:);
        A = pop((2*n)+1:1:3*n,:);
        J = pop((3*n)+1:1:4*n,:);
        Q = pop((4*n)+1:1:5*n,:);
        I = pop((5*n)+1:1:6*n,:);
        H = pop((6*n)+1:1:7*n,:);
        R = pop((7*n)+1:1:8*n,:);
        D = pop((8*n)+1:1:9*n,:);
        Ir = pop((9*n)+1:1:10*n,:);
       
        

        % Set-up the size for the population

        dpop = zeros(size(pop));
        

        

        %% Run the logistic curve
            y = GeneralisedRichard(maxtime);
        
            for m =1:length(y)
                [pi(m)] = y(m)/sum(para.N);
            end


        
        % Age mixing
        AgeMat = (ageMix*(para.a.*para.h)');
        kage=size(AgeMat); %Check the size

        %% SET UP THE DIFFERENTIAL EQUATIONS

        
            for j= 1:soc
                for k= 1:soc
                    SocInf(j,:) = w(j,k)*(para.tau*J(k,:)).*(S(j,:)./para.N(j,:));
                    size(SocInf)

           
               
                
                    
                    dpop(1:1:n,:) = -  SocInf*AgeMat+para.epsilon * R;
                    dpop(n+1:1:2*n,:) =   SocInf*AgeMat*pi(m).*E - (1-pi(m))*para.theta*E - (1-pi(m))*para.rho*E;
                    dpop((2*n)+1:1:3*n,:) = (1-pi(m))*para.rho*E - para.gamma*A;
                    dpop((3*n)+1:1:4*n,:) = (1-pi(m))*para.theta*E - para.delta.*J - para.gamma*J;
                    dpop((4*n)+1:1:5*n,:) = pi(m)*E - 1/para.PHI*para.p*Q - 1/para.PHI*(1-para.p)*para.gamma*Q;
                    dpop((5*n)+1:1:6*n,:) = 1/para.PHI*para.p*Q - para.gamma .*I -para.gamma*I;
                    dpop((6*n)+1:1:7*n,:) = para.gamma .*J+ para.gamma .*I - para.gamma.*H - para.gamma_2*H;
                    dpop((7*n)+1:1:8*n,:) = para.gamma*J + para.gamma*A +para.gamma_2*H +1/para.PHI*(1-para.p)*para.gamma*Q + para.gamma*I - para.epsilon*R;
                    dpop((8*n)+1:1:9*n,:)=  para.gamma.*H;
                    dpop((9*n)+1:1:10*n,:) = 1/para.PHI*para.p*Q;
                end
            end
        
        dpop = dpop(:);

                      
    end
end

The initial conditions are set up like this

ICs = struct('S',S,'E',E,'A',zeros(soc,n),'J',zeros(soc,n), ...
    'Q',zeros(soc,n),'I',zeros(soc,n),'H',zeros(soc,n),'R',zeros(soc,n), ...
    'D',zeros(soc,n),'Ir',zeros(soc,n));

However my output is giving just the re repetition of the first row for the number of times and in each social groups. What am I doing wrong please.

Answer 1

[Torsten]

You write the same data in the Classes structure for each value of i:

     for i = 1:soc
         Classes.S(:,:,i) = pop(:,1:1:n);
         Classes.E(:,:,i) = pop(:,n+1:1:2*n);
         Classes.A(:,:,i) = pop(:,(2*n)+1:1:3*n);
         Classes.J(:,:,i) = pop(:,(3*n)+1:1:4*n);
         Classes.Q(:,:,i) = pop(:,(4*n)+1:1:5*n);
         Classes.I(:,:,i) = pop(:,(5*n)+1:1:6*n);
         Classes.H(:,:,i) = pop(:,(6*n)+1:1:7*n);
         Classes.R(:,:,i) = pop(:,(7*n)+1:1:8*n);
         Classes.D(:,:,i) = pop(:,(8*n)+1:1:9*n);
         Classes.Ir(:,:,i) = pop(:,(9*n)+1:1:10*n);
     end