
dx = 1; %each depth step

Nx = 20; 
Mx = 5; 
Bx = 7; 
Nt = 100; 
totalx = Nx+Bx+Mx; 
bdry1 = totalx - 12; 
bdry2 = totalx - 5; 
dt = (8*60*60)/Nt 

De = 5.*(10.^-10); %molecular diffusivity in epidermis

Cb = zeros(totalx+1,Nt+1); %boundary region

time = [0:28800/Nt:28800];


c = 3.7;

    %this is the loop through space, the sub-loop through time calculates the change
    %for each cycle of the space loop
    for j=1:totalx+1,
    	if (j > bdry1 & j < bdry2) %at defined boundary layer
    		disp ("j = "), disp (j);
    		Cb(j,:) = c.*(sin (0.0001.*time)).*(e.^(-0.000363.*time));
    		Cblast = Cb; %Save the last guess
    		Cb(:,1) = Cblast(:,end); %Initialize the value at t=0 to the last value
    		for i=2:Nt+1, %this is the loop through time, space = constant
    			depth_2D = (Cb(j:end-2,i-1)-2*Cb(j+1:end-1,i-1)+Cb(j+2:end,i-1))/dx^2;
    			time_1D = De*depth_2D;
    			Cb(j+1:end-1,i) = time_1D*dt + Cb(j+1:end-1,i-1);
    			Cb(end,i) = Cb(end-1,i); % Enforce end BC
    		end
    	endif
    end


%Plot the solution





figure(2)
imagesc([0 12],[0 30],Cb); 
title('Top view plot (imagesc)')
colorbar


figure(3)
depth = [0:dx:totalx*dx];
%mesh(time,depth,Cb);
meshc(time,depth,Cb);
title('3D (mesh or meshc)')
colorbar

