function [] = prob_3_20
% 
% Written by:
% -- 
% John L. Weatherwax                2005-04-13
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

beta = 0; 
gamma = 20; 
lambda = 0.5; 
theta_c = 1; 
B = 6; 
D_a = 0.05;
D_a = 0.053;

solinit = bvpinit(linspace(0,1,10),@guess);  
sol = bvp4c(@f,@bcs,solinit,[],beta,gamma,lambda,theta_c,B,D_a);

figure; 
plot( sol.x, sol.y(1,:), '-bo' ); hold on; grid on;
disp( 'Values requested' ); 
deval( sol, [ 0 ] )

function v = guess(x)
% GUESS - 
%   
%v = [ 1; 1 ]; % <- gets first solution
%v = [ 2; 2 ]; % <- gets second solution
v = [ 4; 5 ]; % <- gets third solution


function [res] = bcs(ya,yb,beta,gamma,lambda,theta_c,B,D_a)
% BCS - 
%      

res = [ ya(1)-(1-lambda)*yb(1);
        ya(2)-(1-lambda)*yb(2);
      ];


function [ ode ] = f(x,y,beta,gamma,lambda,theta_c,B,D_a)
% F - 
% y(1)=y(x); y(2)=theta(x);
% 
  
tmp = exp( gamma*y(2)/(gamma+y(2)) );   
  
ode = [  
    D_a*(1-y(1))*tmp;
    B*D_a*(1-y(1))*tmp - beta*(y(2)-theta_c);
      ];