function lv2_dist_3d(r1,k1,n01,r2ini,r2end,r2step,k2,n02,alfa12,alfa21,frequ,intini,intend,intstep,nostep)
%Lotka Volterra system with disturbance
%This function plots response surface of one species to disturbance intensity 
%and growth rate of one species

%r1 - intrinsic growth rates for species 1 
%k1, k2 - carrying capacities
%n01, n02 - initial population sizes
%alfa12, alfa21 - competition coefficients
%frequ - how often to disturb (no of time steps between disturbance events)
%r2ini,r2end,r2step - minimum, maximum growth rate of species 2 and difference between adjacent values
%intini,intend,intstep - minimum, maximum proportion of population of each species 
%   left after disturbance and difference between adjacent values
%nostep - number of steps for simulation



if r2ini >= r2end | r2step<=0 | intini >= intend | intstep<=0
   return;
end

start=min(100,nostep/2);
x=[];
z=[];
for r=r2ini:r2step:r2end
   y=[];
   z1=[];
   for intens=intini:intstep:intend
      x1=lv2_dist(r1,k1,n01,r,k2,n02,alfa12,alfa21,frequ,intens,intens,nostep);
      a=mean(x1(2,start:nostep)/(x1(1,start:nostep)+x1(2,start:nostep)));
      y=[y 1-intens];
      z1=[z1 a];
   end
   x=[x r./r1];
   z=[z;z1];
end
%plot3(x,y,z);
surf(x,y,transpose(z));
xlabel('Relative growth rate of species 2');
ylabel('Disturbance intensity');
zlabel('Proportion of species 2');

return;



