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);
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);
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',[]);
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);
function dpop = diff_socialmodel(t,pop,para)
ageMix = Mix_H + Mix_W + Mix_S + Mix_O;
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,:);
y = GeneralisedRichard(maxtime);
[pi(m)] = y(m)/sum(para.N);
AgeMat = (ageMix*(para.a.*para.h)');
SocInf(j,:) = w(j,k)*(para.tau*J(k,:)).*(S(j,:)./para.N(j,:));
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;