%function [ M ] = implicit_diffusion(x0,b,lambda)

% 

% x0 is the vector of intial postions (at t=0) 

% b is a 2-column vector.  the first column is the left side, the second

% column is the right side

%

 

%find the size of the boundary condition arrays that are passed in.

% we use max(size) to get the length, since size gives both the row and

% column numbers

Length = 1;

Time   = 5;

dt     = 0.01;

dx     = 0.2;

k      = 0.25;

lambda = dt/dx^2;

 

 

%D = 1;

nt = Time/dt ;

nx = Length/dx + 1;

 

%xin = linspace(1,dx*(nx-1),nx)';

%x0 = xin; 

 

x0 = [1 0 0 0 0 0];

 

for i=1:nt

    b(1,i) = 1;

end

 

b(1,1) = 1;

 

 

%nt = max(size(b));

%nx = max(size(x0)); 

 

%allocate the matrix of the solutions

M = zeros(nx,nt);

%allocate the solution matrix for point j+1

A = zeros(nx,nx);

B = zeros(nx,1);

 

%set the initial conditions

M(:,1)  = x0;

k      = 0.25;

 

for j=1:nt-1

 

    %set the x=0 boundary conditions equations at t=j+1

    A(1,1) = 1;

    B(1)   = b(1,j+1);

 

    %set the x=n boundary conditions equation at t=j+1

     A(nx,nx) = 1;

     B(nx)  = M(nx-1,j);

 

    % set the matrix at the interior point nodes at t=j+1

    

    for i=2:nx-1

        A(i,i-1) = -lambda;       

        A(i,i)   = 1+2*lambda+k*lambda*dx^2;

        A(i,i+1) = -lambda;          

        B(i)     = M(i,j);

        %A(1,1) = 1;

        %A(nx,nx-1) = -2*lambda;

        %A(nx,nx) = 1+2*lambda+k*lambda;

    end    

     

    M(:,j+1) = A\B;    

end

M';

 

%M = CrankNicholson(x0,b,lambda);

[t,x] = meshgrid(0:dt:dt*(nt-1),0:dx:dx*(nx-1));

 

surf(t,x,M)

shading interp

 

xlabel('time')

ylabel('distance')