 %Spruce budworm model on [0,L]

 %  u_t = D u_xx + f(u)

 %Neumann BC on [0,L]

function pdexbudworm
m = 0;
L = 1;
x = linspace(0,L,40);
t = linspace(0,0.5,40);
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
%Extract the first solution component as u.
u = sol(:,:,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Surface plot of u(x,t)
figure(1)
surf(x,t,u) 
title('Numerical solution: u(x,t)')
xlabel('x')
ylabel('t')
%Solution profile u(x,*) at last time.
figure(2)
plot(x,u(end,:))
title('u(x,t) at t = xxxx')
xlabel('x')
ylabel('u(x,1) fix this')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %        Defining the problem (m=0)

%  c(x,t,u,u_x) u_t = (f(x,t,u,u_x))_x + s(x,t,u,u_x)

function [c,f,s] = pdex1pde(x,t,u,DuDx)
c = 1;
DD= 0.2;
f = DD*DuDx;

r = 0.5;
q = 9;
s = r*u.*(1-u/q)-(u.^2)./(1+u.^2);
%       u(x,0)=u0(x)

function u0 = pdex1ic(x)
u0=sin(pi*x).^2;


%      p + q*f = 0          x=xl,xr

% where f is the flux f declared in the PDE above

function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
pl = 0;
ql = 1;
pr = 0;
qr = 1;