% Implicit method for u_t = k u_xx

clear;  %  gets rid of old data

L = 10;  % domain size
dt = .1; % time step size
M = 100; % number of time steps

N = 50; % number of spatial mesh points
dx = L/N;

s = dt/dx^2;  %derived parameter that appears in the numerical method

%setup grid

x = [dx: dx: L-dx]';   % x is the column vector x_1, ..., x_N-1
                       % the ' denotes a matrix transpose

% setup initial data:  u is a N-1 size column vector which
% will always contain the solution at the present time step

u = x.*(10-x);  % initial data u = x*(10-x).  The dot before
                % the * indicates to multiply element by element
		% instead of a matrix-matrix product

% A is a N-1 x N-1 matrix which allows the method to be written
% u_m = A*u_m+1.
% the command "diag(vector,k)" constructs a matrix with vector
% on the kth diagonal (k=0 is the central one)

A = (1+2*s)*diag(ones(N-1,1),0) - s*diag(ones(N-2,1),1) - s*diag(ones(N-2,1),-1);

time=0;
for j = 1:M,
   % these commands will plot the solution and pause
   % briefly before proceeding to the next time-step.
   % to make things go slower, increase the pause time

   plot(x,u);
   axis([0 10 -10 25]);
   text(0.1,20,['time=' num2str(time)]);
   pause(.05);

   unext = A\u;        % Solve the equation A*unext =u;
   u = unext;          % now u is equal to the updated solution
   time = time + dt;
end;
