Daisy World Simulation Help?

Hello! I'm doing a simple (or least I thought it would be) Daisy World simulation. However, the area of the black daisies and white daises come out to be the same after 200 time steps, which is definitely wrong. Could anyone look over this code and tell me where I went wrong?

function dy= DaisyWorld(t,y)

%Initial Conditions WhiteArea=0; BlackArea=0; UncoveredArea=1; tspan=[0,200];

%Influences blackAlbedo= 0.25; %albedo measures the reflectivity of a surface

uncoveredAlbedo= 0.5;

whiteAlbedo= 0.75;

deathRate= 0.3;

heatAbsorpFact= 20; %this controls how the local temps of the daisies differ from the average planetary temp

sbConstant= 5.669e-8;

solarFluxConstant= 917;

i=0; for t=1:200; %time steps i=i+1;

end

solarLuminosity= 0.6+(i*(1.2/200)); %solar luminosity increases with time

planetaryAlbedo=(UncoveredArea*uncoveredAlbedo)+(BlackArea*blackAlbedo)+(WhiteArea*whiteAlbedo);

avgPlanetTemp=((solarLuminosity*solarFluxConstant*(1-planetaryAlbedo)/sbConstant)^0.25)-273;

tempBlackLand=heatAbsorpFact*(planetaryAlbedo-blackAlbedo)+avgPlanetTemp;

tempWhiteLand=heatAbsorpFact*(planetaryAlbedo-whiteAlbedo)+avgPlanetTemp;

tempDeadPlanet=((solarLuminosity*solarFluxConstant*(1-0.5)/sbConstant)^0.25)-273;

blackGrowthFact=1-.003265*((22.5-tempBlackLand)^2);

whiteGrowthFact=1-.003265*((22.5-tempWhiteLand)^2);

db= BlackArea*(UncoveredArea*blackGrowthFact-deathRate)+0.001;

dw= WhiteArea*(UncoveredArea*whiteGrowthFact-deathRate)+0.001;

dy=[db;dw]

Run code:

%initial conditions WhiteArea=0; BlackArea=0; UncoveredArea=1; IC=[0,0]; tspan=[0,200];

[t,dy]=ode45(@DaisyWorld,tspan,IC);

plot(t,dy(:,1))%blackarea

legend('blackarea')

hold on;

plot(t,dy(:,2))%white area

legend('whitearea')

title('Area vs. Time')

xlabel('Time')

ylabel('Area')