function budworm clear all; close all; x = [1:0.01:50]; rx = 2*x.^3./(1+x.^2).^2; kx = 2*x.^3./(x.^2-1); figure(1); plot(kx,rx,'linewidth',2); set(gca,'fontsize',20); xlabel('k(x)'); ylabel('r(x)'); grid on; xlim([0 110]); pause; %%%%%%%%% flow plot %%%%%%%%% r = 0.3; k = 30; x = [0:0.01:10]; N = length(x); for i = 1:N dx(i) = r*x(i)*(1-x(i)/k)-x(i)^2/(1+x(i)^2); end; r1 = 0.6; for i = 1:N dx1(i) = r1*x(i)*(1-x(i)/k)-x(i)^2/(1+x(i)^2); end; r2 = 0.1; for i = 1:N dx2(i) = r2*x(i)*(1-x(i)/k)-x(i)^2/(1+x(i)^2); end; figure(2); hold on; h=plot(x,dx,'b',x,dx1,'r',x,dx2,'c'); set(h,'linewidth',2); legend('bistable','outbreak','refuge'); ylim([-0.7 0.7]); set(gca,'fontsize',20); title('k=20'); xlabel('x'); ylabel('dx'); grid; pause; figure(3); hold on; h=plot(x,dx,'b',x,dx1,'r',x,dx2,'c'); set(h,'linewidth',2); legend('bistable','outbreak','refuge'); xlim([0 0.6]); ylim([-0.1 0.1]); set(gca,'fontsize',20); title('k=20'); xlabel('x'); ylabel('dx'); grid; pause; %%%%%%%%% solution %%%%%%%%%%%%%% options = odeset('reltol',1e-6,'abstol',1e-6); N0 = 10; te = 50; par = [k r]; sol = ode15s(@rhs,[0 te],N0,options,par); t1 = sol.x; N1 = sol.y; par = [k r1]; sol = ode45(@rhs,[0 te],N0,options,par); t2 = sol.x; N2 = sol.y; par = [k r2]; sol = ode45(@rhs,[0 te],N0,options,par); t3 = sol.x; N3 = sol.y; figure(3);clf; h=plot(t1,N1,t2,N2,'r',t3,N3,'c'); set(h,'linewidth',2); set(gca,'fontsize',20); xlabel('time'); ylabel('budworms'); legend('bistable','outbreak','refuge',2); grid on; function dN = rhs(x,y,par) N = y(1); k = par(1); r = par(2); dN = r*N*(1-N/k)-N^2/(1+N^2);