%  This function solves the second part of homework problem [H6.1]
%  ---------------------------------------------------------------
%
echo on
diary nsds_hw6_1b.dry
%
%  Get boundary values
%
tmax = 72;
tvec = [  0;  6; 12; 18; 24; 36; 48; 60; 66; 72];
yvec = [  6; 16; 28; 21; 18; 34; 18; 25; 15;  4];
tvec_new = [0:6:72]';
yvec_new = spline(tvec,yvec,tvec_new);
tvec = tvec_new;
yvec = yvec_new;
tvec_new = [0:72]';
yvec_new = spline(tvec,yvec,tvec_new);
%
%  We use the standard Matlab stiff system solver
%
[tvec,yvec] = ode15s('soil');
tvec = tvec';
yvec = yvec';
[n,m] = size(tvec);
%
%  Plot the results
%
plot(tvec_new,yvec_new,tvec,yvec(50,:))
grid on
title('ttt')
xlabel('xxx')
ylabel('zzz')
print -deps nsds_hw6_1b.eps
clf
%
xmax = 50;
dx = xmax/50;;
xvec = [0:dx:xmax]';
tlin = linspace(min(tvec),max(tvec),73);
xlin = linspace(min(xvec),max(xvec),51);
[T,X] = meshgrid(tlin,xlin);
U = [ yvec_new' ; yvec ];
mesh(T,X,U)
title('ttt')
xlabel('xxx')
ylabel('yyy')
zlabel('zzz')
print -deps nsds_hw6_1c.eps
%
diary off
echo off
%
return