%Crank-Nicolson Soln of 1-D Transport Eqn
%Note that Rcell must be less than 2 to avoid
%sawteeth--RSling 3/95;3/00
%jmax is the max number of nodes (including Boundaries)

% Plot the data first:
load Elev.dat;
x=(0:2.5:49);
plot(x,Elev, '*');
hold;
axis([0 50 1000 2200]);

pause;
%Set some constants

jmax=;
tmax=;
delx=;
delt=;
u=;

%Zero out, define the length of some arrays, and
%change some vectors from row to column

for j=1:jmax;
k(j)=;
qw(j)=;
dkqwdx(j)=;
C(j)=dkqwdx(j)*delt/delx;
D=k(j)*qw(j);
s(j)=D*delt/delx^2;
Rcell=C(j)*delx^2/(delt*D)

%Initial Condition
T(j)=2077-19.0*j;
end
for j=1:jmax-2;
TT(j)=0.0;
B(j)=0.0;
end
B=B';
T=T';

%Fill in coefficient matrix

for j=1:jmax-2;
for i=1:jmax-2;
if i==j,
     A(i,j) = 2.0*(1+s(j));
   elseif (i-j) == 1,
     A(i,j)=-(s(j)+0.5*C(j));
   elseif (i-j) == -1,
     A(i,j)=-(s(j)-0.5*C(j));
   else
     A(i,j)=0;
end
end
end

%Start time loop

for k=1:tmax;

% Fill in RHS

 for j = 2:jmax-1;
     B(j-1)=(s(j)+0.5*C(j))*T(j-1)+2.0*(1-s(j))*T(j);
     B(j-1)=B(j-1)+(s(j)-0.5*C(j))*T(j+1)+u*delt;
 end
 B(1)=B(1)+(s(1)+0.5*C(1))*T(1)+u*delt;
 B(jmax-2)=B(jmax-2)+(s(jmax)-0.5*C(jmax))*T(jmax)+u*delt;

%Now comes the happy moment: invert matrix

TT=A\B;

% Assign new values of T to appropriate nodes
 
 for j=2:jmax-1;
   T(j)=TT(j-1);
 end
 
%Set BC's

 T(1)=T(2);
 T(jmax)=1090.0;
 
%Let's look at it



   plot(x,Elev, '*');
   axis([0 50 1000 2200]);
   plot(T(1:jmax));
   hold;
   drawnow;
 end
