%  This code solves the heat equation 
%
%      u_t = u_xx ,  0 <= x <= 1,  t > 0,
%
%  with IC:  u(x,0) = x(1-x),  0 <= x <= 1,
%  and BC:   u(0,t) = u(1,t) = 0,  t > 0.
%
%  User enters:  n = no. of interior mesh points in x,
%               dt = time step,
%             tmax = time at which solution is desired.
%
%  For the homework exercise, set n=19 (so that dx=1/(n+1)=.05) and
%  tmax = 0.1.  Try  dt = .25*(.05)^2,  .5*(.05)^2,  and (.05)^2.
%  The first two should give reasonable answers; the third dt value
%  should be unstable and give wild oscillations.
%

n = input(' Enter n: ');
dt = input(' Enter dt: ');
tmax = input(' Enter tmax: ');

%  Set initial values.

dx = 1/(n+1);
u = zeros(n,1);
for i=1:n,
  x = i*dx;
  u(i) = x*(1-x);
end;
uold = u;

%  Loop over time steps.

for m=1:tmax/dt,
  t = m*dt;
  u(1) = uold(1) + (dt/dx^2)*(uold(2)-2*uold(1));
  for i=2:n-1,
    u(i) = uold(i) + (dt/dx^2)*(uold(i+1)-2*uold(i)+uold(i-1));
  end;
  u(n) = uold(n) + (dt/dx^2)*(-2*uold(n)+uold(n-1));
  uold = u;
end;

%  Plot solution.

plot(dx*[1:n]',u)