function [P,P_ratio] = run_hypergraph_mean_field_approximation(i_0,beta,delta,I,T,dt,K,fun)
% i_0: initial proportion of infectious individuals
% beta, delta: rate of infectionm resp. recovery
% I: nxm incidence matrix, where n is the number of nodes and m the number
% of hyperedges
% f: non-linear dependence on the number of infectious neighbours
% T: finite time
[n,m]=size(I);
n_steps=floor(T/dt);
P=zeros(n,n_steps+1);
P_ratio=zeros(1,n_steps);
P(:,1)=i_0;

for t=1:n_steps    
    %for i=1:n
            P_t=P(:,t);
            %P_t(i)=0;
            P(:,t+1) = P(:,t)+dt*(beta*sum(I.*f(sum(P(:,t).*I,1),sum(I,1),K,fun),2).*(ones(n,1)-P(:,t))-delta*P(:,t));
            %P(i,t+1) = P(i,t)+dt*(beta*sum(I(i,:).*f(sum(P_t.*I,1),sum(I,1),K,fun),2)*(1-P(i,t))-delta*P(i,t));
    %end
            P_ratio(t)=sum(P(:,t+1))/sum(P(:,t));
end
end