function [out1,out2,out3] = soil(t,x,flag)
%
%  State-space model of [H6.1]
%
%  First, get boundary value at current time
%
tvec = [  0;  6; 12; 18; 24; 36; 48; 60; 66; 72];
uvec = [  6; 16; 28; 21; 18; 34; 18; 25; 15;  4];
tvec_new = [0:6:72]';
uvec_new = spline(tvec,uvec,tvec_new);
tvec = tvec_new;
uvec = uvec_new;
tvec_new = [0:72]';
u0 = spline(tvec,uvec,t);
%
%  Next, compute the diffusion coefficient
%
lambda = 0.004;
rho = 1.335;
c = 0.2;
dx = 1;
sigma = 3600*lambda/(rho*c*dx*dx);
%
%  Get the matrices of the linear state-space model
%
A = sigma*(diag(-2*ones(50,1)) + diag(ones(49,1),1) + diag(ones(49,1),-1));
A(50,49) = sigma*2;
b = [ sigma*u0 ; zeros(49,1) ];
%
if nargin < 3 | isempty(flag)
  out1 = A*x + b;
 else
  switch(flag)
   case('init')
    tmax = 72;
    smax = 1;
    ds = smax/50;
    out1 = [0:tmax];
    out2 = 6*ones(50,1);
    out3 = odeset('RelTol',1.0e-5,'MaxOrder',4,'NormControl','on')
   case('jacobian')
    out1 = A;
   otherwise
    error(['Unknown flag ''' flag '''.']);
  end
end
%
return