% Fishery Management Plan
% Modify line 42 with your own function h(t).  
% If your harvesting is dramatically too high, you'll get an error like "Failure at t=..." 

T = [0,25];
U = chebop(T);
U.numVars = 1;
U.op = @fish_differential_equation;
U.lbc = 2500;
y = U\0;
y = max(y,0);
h = chebfun(@harvest_rate,T).*(y>0.5);

sus = 2000 * (floor(y(T(2))) >= 1500);
sal = sum(14*chebfun('exp(-0.1*t)',T).*h);
J = sus + sal;

fprintf('\n Total fish eaten: %0.0f\n',sum(h))
fprintf(' Final fish population: %0.0f\n',floor(y(T(2))))
fprintf(' Sustainability bonus: %0.2f\n Harvested fish value: %0.2f\n total J: %0.2f\n\n',sus,sal,J)

figure(1)
clf
plot(y,'LineWidth',2)
ylim([0,10000])
ylabel('Fish population')
xlabel('Time in years')
title('Fishery population forecast')

figure(2)
clf
plot(h,'Linewidth',2)
ylabel('Harvest rate (fish per year)')
xlabel('Time in years')
title('Fishery harvesting plan')

function dydt = fish_differential_equation(t,y)
    dydt = -diff(y) + (0.2*y.*(1-y/10000) - harvest_rate(t));
end

function h = harvest_rate(t)
    h =  100 + 10*t ;
end