Cant get my output to display both black and white daisies?

Hi,

I am having trouble with a model run. I am testing the daisyworld hypothesis and I have been able to show how the black and white daisies both grow seperately, but I cannot get my graph to show how they grow simultaneously. Can anyone help?

Below is my model to show the black daisy growth, to get the white daisies I just change the black albedo and IC to 0. I have tried everything I can think of to get both showing.

function test
IC = [0.3 0]; %initial daisy coverage percentages, black & white respectively  
tspan = [0 3000]; %solve for values of t  
L = @(t) t/range(tspan)*2 + 0.4; %luminosity function (has to stay between 0 and 2)   
[T,Y] = ode45(@(t,y)daisyworld(t,y,L),tspan,IC); % solves equation, need capital T and Y   
Lum= L(T); % luminosity at time t  
PT = Temp(Lum,Y); % planetary temperature over time  
subplot(211), plot(T,Y), xlabel('T, Time'), legend('Black daisies','White daisies')  
subplot(212), plot(Lum,PT)
xlabel('L(T), Luminosity'), ylabel('PT, Plantary Temperature')

function PT = Temp(Lum,y)  
if isvector(y), y = y(:)'; end % create y row vector  
solar = 917; % solar constant  
SB = 5.67e-8; % Stefan-Boltzmann constant;  
Ab = 0.25; % Albedo of black daisies  
Aw = 0; % Albedo of white daisies  
Ag = 0.5; % Albedo of bare ground  
A = ((1-y(:,1)-y(:,2))*Ag)+(y(:,1)*Ab)+(y(:,2)*Aw);  % albedo  
PT = ((((solar*Lum)/SB).*(1-A)).^(1/4))-273; % plantary temperature 

function dydt = daisyworld(t,y,L)  
dydt = [0;0];  
Lum = L(t); % calculate Luminosity at time t  
PT = Temp(Lum,y);  
ab = y(1); % area colonized by black daisies  
aw = y(2); % area colonized by white daisies  
B = 1-0.003265*((22.5-PT).^2); % growth rate (beta value)  
g = 0.3; % death rate  
dydt(1) = ab*(((1-ab-aw)*B)-g); % black daisy formula  
dydt(2) = aw*(((1-ab-aw)*B)-g); % white daisy formula

Any help would be HUGELY appreciated