function [dxdt]=fishery(t,x)
%Simulation of a fishery
global a b c r K k ki kt xr
r=0.2;K=100;a=0.1;b=1;c=1;
k=0.5;ki=0.5;kt=0.5;t1=20;
f=r*x(1)*(1-x(1)/K);
xr0=c/a/b;
xr=(K+c/a/b)/2;
v=k*(xr-x(1))+x(2);
u=max(0,-v);
h=a*x(1)*u;
g=b*h-c*u;
dxdt=[f-h;ki*(xr-x(1))+kt*(u-v)];
