%
%   HW #1, problem 1.10.2
%
% Solve an ordinary differential equation using Euler's method
%
%   Jon Collis
%   15 January 2009
%

% Just make sure we don't have extra stuff confusing things
clear;


% Define some constant parameters of the problem
Gamma = 4*10^(-3);
mu = 10^(-3);
cstar = 7.52*10^(-7);
c0 = 1.05*cstar;
xstar = .05;
k = 5*10^7;

c1 = c0 + mu*xstar^3;

% The book gives these values, so we verify them here
k*mu
k*cstar
k*c1


% Start doing some Euler method stuff.

tmax = 0.5;  % maximum time in seconds

% initial values
t0 = 0;
y0 = xstar;

t(1) = t0; % gonna save time and y values here, though don't have to
y(1) = y0;

h = 0.0001;  % This is the grid spacing



imax = tmax/h;

for j=1:imax
    j
    t(j+1) = t(j) + h;
    y(j+1) = y(j) + h*(k*(c1 - mu*y(j)^3 - cstar*exp(Gamma/y(j))));
end


'Our solution'
figure(14); clf;
plot(t,y,'LineWidth',2);
xlabel('time (s)');
ylabel('crystal size (\mum)');
title('1.10.2')
set(gca,'XLim',[0 0.5]);
set(gca,'YLim',[0.045 0.051]);


% Just for fun, draw the xi_2 line
hold on
plot([0 0.5],[.0453479 .0453479],'r-')
hold off

% let's verify that the solution is greater than xi_1 and xi_2
% from the book, xi_1 and xi_2 are about
% .0216736 and .0453479
%
max(y)
min(y)  % Note that the minimum value for this case is going to be greater than both xi_1 and xi_2