clear

%% Parameters
r1 = 5; % gallons/minute (rate in and out during Stage 1) 
r2 = 5; % gallons/minute (rate in and out during Stage 2) 
c0 = 1.3e-5; % ounces per gallon of water (initial concentration)
c1 = 1e-7; % ounces per gallon of water (concentration flowing in during Stage 1)
c2 = 5e-5; % ounces per gallon of water (concentration pumped in during Stage 2)
V = 84000; % gallons (volume of pool)
T = 60*24; % minutes (duration of Stage 1)

%% Stage 1: Right-hand side
RHS1 = @(Q,t) c1*r1 - Q*r1/V;

%% Stage 2: Right-hand side
RHS2 = @(Q,t) c2*r2 - Q*r2/V;


%% Solve the differential equation of Stage 1
tspan1=[0,T];
h=0.1; % Step size
Q1ini = c0*V; % Initial condition for Stage 1
[Q1Numerical,t1] = Euler(RHS1,tspan1,Q1ini,h);

%% Solve the differential equation of Stage 2
tspan2=[T,2*T];
Q2ini = interp1(t1,Q1Numerical,T); % Initial condition for Stage 2
[Q2Numerical,t2] = Euler(RHS2,tspan2,Q2ini,h);

%% Exact solutions
Q1Exact=V*(c0-c1)*exp(-r1*t1/V)+V*c1; % Stage 1
Q2Exact=((c0-c1)*exp(-r1*T/V)+(c1-c2))*V*exp(-r1*(t2-T)/V)+V*c2; % Stage 2

%% Plots
figure(1)
plot([t1;t2],[Q1Numerical-Q1Exact;Q2Numerical-Q2Exact])
hold on
xlim([0 2*T])
%% Add title 
title('Numerical error Q_{numerical}-Q_{exact}')
%% Add legend 
legendHandle=findobj(gcf, 'Type', 'Legend');
if length(legendHandle)==0; 
    legendText{1}=strcat('h=',num2str(h));
else
    legendText=legendHandle.String;
    legendText{end}=strcat('h=',num2str(h));
end
legend(legendText,'location','southeast')
%% Add lines
if h==0.01, 
    yline(0);
    xline(T);
end
legend(legendText,'location','southeast')
%% Add axis labels
xlabel('minutes')
ylabel('ounces')


%% Table of errors for various step sizes
for h = 0.1:-0.01:0.01
    %% Numerical solutions
    [Q1Numerical,t1] = Euler(RHS1,tspan1,Q1ini,h); % Stage 1 numerical solution
    Q2ini = interp1(t1,Q1Numerical,T); % Initial condition for Stage 2
    [Q2Numerical,t2] = Euler(RHS2,tspan2,Q2ini,h);
    %% Exact solutions
    Q1Exact=V*(c0-c1)*exp(-r1*t1/V)+V*c1; % Stage 1
    Q2Exact=((c0-c1)*exp(-r1*T/V)+(c1-c2))*V*exp(-r1*(t2-T)/V)+V*c2; % Stage 2
    %% Errors
    format short e
    disp([h,max([Q1Numerical-Q1Exact;Q2Numerical-Q2Exact])])
end

%% Auxiliary function
function [y,t] = Euler(f,tspan,y0,h)
%EULER Solves differential equations using Euler's method. 
% 
%   The differential equation to be solved is of the form
%   y'(t) = f(y,t),
%   subject to the initial condition
%   y(0) = y0.
%
%   The command 
%     [y,t] = Euler(f,tspan,y0,h)
%   solves the problem over the time interval tspan (a two-element vector) 
%   with the step size h. 

N=ceil((tspan(end)-tspan(1))/h); % Number of grid points
y=zeros(N,1); % Predefine y
y(1)=y0; % Apply the initial condition
for i=1:N
    y(i+1) = y(i) + h*f(y(i),i*h); % Integrate in time using Euler's method
end
t=tspan(1)+h*(0:1:N);
t=t(:);

end