B = 0.0065;
lambda = 0.08;
l = .0001;
p = input('enter p: ');
tstep = .000001;
i = 2;
P(1) = 1;
C(1) = lambda./B;
tval(1) = 0;

for t=0:tstep:.001
    dP = ((p-B)./l) * P(i-1) + lambda * C(i-1);
    dC = (B./l) * P(i-1) - lambda * C(i-1);

    if P(i-1) + dP * tstep>5, break, end
    if P(i-1) + dP * tstep<.000001, break, end

    P(i) = P(i-1) + dP * tstep;
    C(i) = C(i-1) + dC * tstep;

    tval(i)=t;
    i=i+1;
end

w1 = 1./2 .* ( - (lambda + B./l - p./l) + ( (lambda + B./l - p./l).^2 + (4.*p.*lambda)./l ).^(1/2) );
w2 = 1./2 .* ( - (lambda + B./l - p./l) - ( (lambda + B./l - p./l).^2 + (4.*p.*lambda)./l ).^(1/2) );
A1 = (p./l - w2) ./ (w1 - w2);
A2 = (w1 - p./l) ./ (w1 - w2);
B1 = A1 ./ (w1 + lambda);
B2 = A2 ./(w2 + lambda);

x = 0:0.0001:t;
P2 = A1 .* exp(w1 .* x) + A2 .* exp(w2 .* x);
C2 = B1 .* exp(w1 .* x) + B2 .* exp(w2 .* x);

plot(tval, P, 'b', tval, C, 'r', x, P2, 'c', x, C2, 'k')
axis([0 t * 1.1 0 13])
legend('P(t)', 'C(t)', 'P(t) exact', 'C(t) exact', 'location', 'EastOutside')
title(strcat('Point Reactor Kinetics Equations with One Delayed Neutron for p= ', num2str(p), ' $'))
xlabel('time')
ylabel('P(t) and C(t)')