function P3()
% Global Warming Simulation
global fCount
close all
% The CO2 Source function...
year = [1000, 1850, 1950, 1980, 2000, 2050, 2080, 2100, 2120, 2150, 2225,...
        2300, 2500, 5000]';
rate = [0, 0, 1.0000, 4.0000, 5.0000, 8.0000, 10.0000, 10.5000, ...
        10.0000, 8.0000, 3.5000, 2.0000, 0, 0]';
f=@(t)pchip(year, rate, t);

% Initial condition and time interval of interest...
tspan=[1000 5000];
y0 = [1.00; 2.01; 2.23; 2.20; 2.26];
opts=odeset();

% ode23 solution with default tolerances...

fCount = 0;
[tOut, yOut] = ode23(@derivative, tspan, y0, opts, f);
rate =f(tOut);
figure
subplot(2,1,1);
plot(tOut,rate,'k');
%legend('Fossil Fuel Source Term');
ylabel('fuel');
title('Carbon in the Atmosphere and Ocean','fontsize',14);
subplot(2,1,2);
title(sprintf('Method = ode23    f-evals = %1d',fCount),'fontsize',14)
hold on;
plot(tOut, yOut(:,1),tOut, yOut(:,2),tOut, yOut(:,3),'LineWidth',2);
%legend('atmosphere','shallow','deep');
xlabel('time (yr)');
ylabel('carbon');


% ode23s solution with default tolerances...

fCount = 0;
[tOut, yOut] = ode23s(@derivative, tspan, y0, opts, f);
rate =f(tOut);
figure
subplot(2,1,1);
plot(tOut,rate,'k');
%legend('Fossil Fuel Source Term');
ylabel('fuel');
title('Carbon in the Atmosphere and Ocean','fontsize',14);
subplot(2,1,2);
title(sprintf('Method = ode23s    f-evals = %1d',fCount),'fontsize',14)
hold on;
plot(tOut, yOut(:,1),tOut, yOut(:,2),tOut, yOut(:,3),'LineWidth',2);
%legend('atmosphere','shallow','deep');
xlabel('time (yr)');
ylabel('carbon');


function dydt = derivative(t,y,f)
global fCount
d=8.64;
mu1=4.95*1e2;
mu2=4.95*1e-2;
v_s=0.12;
v_d=1.23;
w=1e-3;
k1=2.19*1e-4;
k2=6.12*1e-5;
k3=0.997148;
k4=6.79*1e-2;

p = y(1);
sigma_s=y(2);
sigma_d=y(3);
alpha_s=y(4);
alpha_d=y(5);

h_s=(sigma_s-(sigma_s^2-k3*alpha_s*(2*sigma_s-alpha_s))^(1/2))/k3;
c_s=(alpha_s-h_s)/2;
p_s=k4*h_s^2/c_s;

dp=(p_s-p)/d+f(t)/mu1;
dsigma_s=1/v_s*((sigma_d-sigma_s)*w-k1-(p_s-p)/d*mu2);
dsigma_d=1/v_d*(k1-(sigma_d-sigma_s)*w);
dalpha_s=1/v_s*((alpha_d-alpha_s)*w-k2);
dalpha_d=1/v_d*(k2-(alpha_d-alpha_s)*w);

dydt=[dp; dsigma_s; dsigma_d; dalpha_s; dalpha_d];
fCount = fCount+1;

% Correct CO2 source plot: 1pt
% Correct plots of p and sigma: 2pt
% Correct f-evals: 2pt