
function salamander


% X0 = initial conditions
% X0 =[F0;B0;Y0;J0;C0;K0;A0;P0];

X0=[1;0.1;0;0;0;0;0.1;0];

% T is the time values
% Xapp is a matrix with all the solutions at the various time 
% steps
% eqns holds the nonlinear system
% [0 1000] is the time interval we care about
% X0 is the intial condition

[T,Xapp] = ode45(@eqns,[0 1000],X0);

figure(1)
plot(T,Xapp(:,1),T,Xapp(:,2),T,Xapp(:,3),T,Xapp(:,4),T,Xapp(:,5),T,Xapp(:,6),T,Xapp(:,7),T,Xapp(:,8))
title('plot')
xlabel('Time')
ylabel('Population')
legend('F','B','Y','J','C','K','A','P')


disp(sprintf('At last time step F = %g',Xapp(end,1)))
disp(sprintf('At last time step B = %g',Xapp(end,2)))
disp(sprintf('At last time step Y = %g',Xapp(end,3)))
disp(sprintf('At last time step J = %g',Xapp(end,4)))
disp(sprintf('At last time step C = %g',Xapp(end,5)))
disp(sprintf('At last time step K = %g',Xapp(end,6)))
disp(sprintf('At last time step A = %g',Xapp(end,7)))
disp(sprintf('At last time step P = %g',Xapp(end,8)))


return

function dy = eqns(t,y)
dy = zeros(8,1);
F = y(1);
B = y(2);
Y = y(3);
J = y(4);
C = y(5);
K = y(6);
A = y(7);
P = y(8);

dy(1) = 10*F*(1-F) - 2*B*F - 2*Y*F;
dy(2) = (1*2*B*F - (.5*P+ .5*A + .1*K + .5*J +.1*C)*B -.001*B);
dy(3) = (.3*(P*(.1*Y+.5*B) + .5*A*B) +...
        .3*(.1*K*B + K*(.1*P+.1*A+.2*J +.1*Y) + .3*.2*K*C)-...
         (.1*K + .5*C + .1*P)*Y - .05*2*Y*F - .001*Y);
dy(4) = .9*.05*2*Y*F - (.5*C + .2*K)*J - .02*.5*J*B -.001*J;
dy(5) = .1*.05*2*Y*F - .2*C*K - .02*C*(.5*Y + .1*B + .5*J) - .001*C;
dy(6) = .02*C*(.5*Y + .1*B + .5*J) - .001*K;
dy(7) = .02*.9*.5*J*B -.1*K*A- .001*A;
dy(8) = .02*.1*.5*J*B - .1*K*P- .001*P;




return

