%  Solves the time-dependent heat equation in one spatial
%  dimension using the forward Euler method.
%  Plots the solution at each timestep and returns solution
%  at t=.1 in vector usave, so that this can be compared
%  with the solution computed using other mesh sizes.
%  (nsteps must be a multiple of 10 in order to get the solution 
%  at t=.1).

%  Set up spatial grid.
n = input(' Enter number of subintervals in space n: ');
h = 1/n;
x = [h:h:1-h]';

%  Set tmin and tmax and timestep size dt.
tmin = 0; tmax = 1;
nsteps = input(' Enter number of time steps nsteps: ');
dt = (tmax-tmin)/nsteps;

%  Initialize uk and plot it.
uk = x.*(1-x);  % This sets uk(i) = x(i)*(1-x(i)) for i=1,...,n-1.
plot(x,uk), pause(1), hold on;

mu = dt/h^2;
A = sparse((1-2*mu)*eye(n-1,n-1));
for i=1:n-2, A(i,i+1)=mu; A(i+1,i)=mu; end;

for k=1:nsteps,
  t = k*dt;
  ukp1 = A*uk;
    if abs(t-.1) < abs(t+dt-.1) & abs(t-.1) < abs(t-dt-.1),  % for t=.1,
      usave = ukp1;                                          % save solution
    end;
  plot(x,ukp1), pause(.1)
  uk = ukp1;
end;
hold off;