function [x, H] = diffuse(dx, Hinitial, k, finaltime, pt, Harb)
%diffuse change topography by diffusive processes
%diffuse(dx, Hinitial, k, finaltime, pt, Harb) (pt and Harb optional)
% Input variables:
%  Scalar
%	dx = distance between nodes
%	k = diffusion constant
%	finaltime = final time of model run
%	pt = number of timesteps between plots
%
%  Vectors
%	Hinitial = initial elevation vector.  The length defines the number of 
%		calculation nodes
%	Harb = vector of arbitrary elevation changes that occur over the 
%		time run
%
% Output variables:
%  Vectors
%	x = xlocations of nodes
%	H = elevations at nodes
%
% Internal variables
%	lambda = stability criterion
%		lambda = dt*k/dx^2,
%		which has to be less than 0.5,
%		so dt = (0.1*dx^2)/k
%	dt = timestep
% 	time = current time of model run
%	i = counter
%	nt = number of time steps( nt = finaltime/dt)

%Start with input
fprintf('Input control variables\n')
fprintf('dx = %.3f k = %.3f\n',dx, k)
dt = (0.1*dx^2)/k;
nt = ceil(finaltime/dt)+1;
fprintf('given those, dt = %.3f, and number of timesteps = %.3f\n',dt, finaltime/dt)

% manage inputs
if nargin < 6, Harb = zeros(1,length(Hinitial));end
if nargin < 5, pt = nt;end

x = [0:dx:((length(Hinitial)-1)*dx)];
H = Hinitial;
fprintf('plotting initial profile and setting up plots\n');
clf;
subplot(2,1,1);
plot(x,H,'k',x,H,'+')
title('Profiles with time');
xlabel('distance along profile,m');
ylabel('elevation,m');
hold on;

subplot(2,1,2);
title('Flux with time');
xlabel('distance along profile,m');
ylabel('flux,m^2');
hold on;

%now erode and deposit
fprintf('running the model and plotting\n');
i = 1;
xhalfs = [x(2)-dx/2:dx:x(length(x))]; 
time = 0;
while time <= finaltime,
%	fprintf('time = %.3f\n',time);
	Q = k.*diff(H)./dx;
		if i==pt
			subplot(2,1,2)
			plot(xhalfs,Q,'k+');	
			drawnow;
		end
	dhint = (diff(Q)./dx).*dt;
	dh = [0 dhint 0];
	H = H + dh + (Harb./(nt+1));
		if i==pt
			subplot(2,1,1);
			plot(x,H,'r');
			drawnow;
%comment the next two lines out if you don't like how the time is being written on the side of the profile plot

			s = sprintf('time = %.3f', time);
			text(x(length(x-2)),max(H),s);
			drawnow;
			i=1;
		end
        time = time +dt;
	i = i+1;
end
%final plot
 subplot(2,1,1);
 plot(x,H,'r');
 s = sprintf('time = %.3f', time-dt);
 text(x(length(x-2)),max(H),s);

hold off;