function [h]=thetarule(theta,heq,dt,T,recoveryTime,hinitial)
% theta         =   theta
% heq           =   equilibrium profile
% dt            =   1
% T             =   halveringstijd oid
% recoveryTime  =   months of recovery
% hinitial      =   initial profile
%
% theta=1;heq=A';dt=1;T=10;recoveryTime=20;hinitial=B'
h(1,:)=hinitial;
for n=1:recoveryTime
    h(n+1,:)=(dt*(heq/T)+((1-((dt*(1-theta))/T))*h(n,:)))/(1+((dt*theta)/T));
end

% for t=1:300
%     plot(h(t,:)), hold on
%     plot(A,'g')
%     pause(0.2);hold off
% end