%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This is in thé main script:
field = [1,0,1 ; 1,1,1];
%in my script i loop through each element of the field-matrix and excute code afterwards but to make it easier to %understand i left that out since its not important for the question
idx1 = 1;
idx2 = 2;
startA = 10;
A = 10;
currbugs = [1;1;1;1;1;1;7];
tspan = [0,1];
y0 = [startA];
[t,y]=ode45(@(t,y)Aph(t,y,field,idx1,idx2,currbugs),tspan,y0);
plot(t,y)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%This is the function which is in a seperate file:
function sol=Aph(t,y,field,idx1,idx2,currbugs)
rgr = 10;
A = y(1);
P = (PRL1(field,idx1,idx2,A)*currbugs(2,1))+(PRad(field,idx1,idx2,A)*currbugs(7,1));
sol = rgr*A*(1-(A/2500))-P;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Within the calculation of the ODE45 equation other function PRL1 and PRad are called but they %are also in seperate files:
function out = PRL1(field,idx1,idx2,A)
sr = 39.1;
ht = 10;
if field(idx1,idx2) == 1
SAI = 1;
else
SAI = 0.6;
end
out = (sr*(A/SAI))/(1+sr*ht*(A/SAI));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%and the second function:
function out = PRad(field,idx1,idx2,A)
sr = 0.0139;
ht = 0.0049;
if field(idx1,idx2) == 1
SAI = 1;
else
SAI = 0.6;
end
out = sr*(A/SAI)/(1+sr*ht*(A/SAI));
end
