% simple Matlab program to calculate rainsplash erosion

C     = 3e-4;         %rainsplash coeff.in years
delx  = 1;            %bins are 1 m wide
delt  = 10;           %time steps are 10 year
total_time = 10000;   % run for 10000 time steps

z = [3 3 3 3 2.25 1.5 .75 0 0 0 0];
                % a list of the current z elevations at x=0 to x=10 {11 nodes}
z_new = z;      % intialize z_new, the place where we store the next time profile

for t=0:delt:total_time;  % the code "0:10:total" means go from t=0 to the total time in delt (10) year steps

    z=z_new;            % update profile to new 'z's, and continue
    for x=2:1:10;       % unfortunately matlab doesn't allow a '0' vector index
                        % so that we have to start from '1' and go to '11', and then we recognize that
                        % the two end nodes are fixed 'BCs', so only calculate interior points
        z_new(x)=z(x) + (C*delt/(delx*delx))*(z(x+1)-2*z(x)+z(x-1));   %this is our equation from class
    end
    z_new(1)= 3;
    x_new(11)=0;              %boundary conditions
    if (mod(t,500) == 0)      % only plot every 500 years
        plot((0:1:10),z);	%again the code '0:1:10' means 0,1,2,3....,10, in shorthand
        hold on
    end
end
axis equal
title('Rainsplash on a fault scarp');
xlabel('meters');
ylabel('Your labels will be diferent');