% samplefd_time_dep.m
% Mousumi Roy, Nov 2014
clear all; 
% finite difference solution of dy/dt = kappa*d^2y/dx^2;  
% As before, lets use the boundary condition y=0.567 at x=0;
% and set y=1 at x=1.
% and lets use an initial condition that is y(x,0) = 0.  Solve 
% for y(x,t) using finite difference (specifically, an explicit 
% method of finite difference, which may or may not be stable):
% So, if i is the space index and j is the time index, then
% y(i,j+1) - y(i,j) =  (dt*kappa/dx)*(y(i+1,j) 
%                     + y(i-1,j)) - (2*dt*kappa/dx)*y(i,j)
% Now rearrange so that we get the unknown quantity y(i,j+1) isolated
% y(i,j+1) = y(i,j)(1 - 2*dt/(dx*dx)) + (dt*kappa/dx^2)*(y(i+1,j)+ y(i-1,j))
% We will turn this into a matrix equation, as we did the case of the 
% steady-state code you wrote for your PS lite.

% lets define a domain and a discretization
dx = 0.1;        
x  = [0:dx:1]'; % domain goes from 0 to 1
n  = length(x);  % there are n points

leftbc  = 0.567;
rightbc = 1;

%define a starting timestep dt
kappa = 1;
dt = 0.00001;
%define a maximum time, maxt
maxt = 3;
maxj = maxt/dt; %maximum number of iterations
nout = 1000; %plot output every nout iterations
% start time
t=0;

% define useful coefficients
C1 = 1 - 2*kappa*dt/(dx*dx); 
C2 = kappa*dt/(dx*dx);

%define a useful matrix that only is dependent on x
X2 = -2*dt*x.*x;
X2(1) = 0;

% first set up initial condition and apply boundary condition at x=0; 
% y(0,1) = 0.567
% y(1,1) = x(1);
% y(2,1) = x(2);
% etc.
yinit = zeros(n,1); 
yinit(1,1)   = leftbc;
yinit(end,1) = rightbc;

% Here are discretized system of equations
% for an arbitrary time j+1, given known quantities at time j:
% y(i,j+1)  =  y(i,j)*C1  +  y(i+1,j)*C2  +  X2;
% or in matrix form for n=5 and j=1:
% [ 1   0   0   0   0 ] [y(1,1)]         = y(1,2)
% [ C2  C1  C2  0   0 ] [y(2,1)]         = y(2,2)
% [ 0   C2  C1  C2  0 ] [y(3,1)]         = y(3,2)
% [ 0   0   C2  C1  C2] [y(4,1)]         = y(4,2)
% [ 0   0   0   0    1] [y(5,1)]         = y(5,2)
% note that the coefficients for the first and last rows have to reflect
% the boundary conditions and either a forward or backward difference for 
% the spatial derivative.
% Let A be the matrix on the left - made of constants and zeros
% let ynew be the matrix on the right
% we get A*y  = ynew

A = C1*diag(ones(n,1)) + C2*diag(ones(n-1,1),-1) + C2*diag(ones(n-1,1),1); 
A(1,1) = 1;
A(1,2) = 0;
A(end, end-1)=0;
A(end,end) = 1;

% now do a time loop to step forward from initial condition
y = yinit;
figure(1); clf 
    plot(x,y(:,1),'r.'); hold on
    plot(x,y(:,1),'k--'); title(['time = ', num2str(t)])

for j = 1:maxj
    t         = t + dt;
    ynew      = A*y ;
    ynew(1,1) = leftbc;
    %some measure of change in y
    del(j,1)  = sqrt( sum( (ynew - y).^2 )/n );
    y         = ynew;
    tt(j,1)   = t;
    %plot output every nout iterations
    if mod(j,nout) == 0
        figure(1); clf
            plot(x,y,'r.'); hold on
            plot(x,y,'k--'); title(['time = ', num2str(t)]); %pause(.5)
    end
end

figure(2); clf
%check steady state
semilogy(tt, del, 'k-'); hold on
