global nu sigma k kt L q h lambda

nu = 0.1;
%nu = 0.04;
sigma = 1.2;
k = 0.4;
kt = 0.4;
L = 10^6;
q = 100;
h = 10;
lambda = q*sigma/h;

[T,Y] = ode23('func_camp_2d',[0,2500],[92.366,2]);  

for i=1:size(Y(:,1))
   bbeta(i) = q*sigma*phi(Y(i,1),Y(i,2))/kt; %construct beta according to
                                             %the quasi-steady state approx
end

figure(1);
hold on;
plot(T,Y(:,1),'b');
plot(T,bbeta,'r');
plot(T,Y(:,2),'g');
xlabel('t');
ylabel('\alpha,\gamma');

