Unverified solution in MATLAB

% Possible solution to Kees' problem in MATLAB
%
% Based on eqns. presented on
% http://websrv.cs.umt.edu/isis/index.php/Kees'_assignment
% Possibly useful as a solution template for those familiar with MATLAB.  
% Note possible problem in lines 55-58.  
 
clear all
close all
clc
 
% Define constants.  
g = 9.8;            % m/ s^ 2; gravitational acceleration
rho = 920;          % kg/ m^ 3; density of ice
delta_x = 100;      % m; distance step
x_max = 25* 10^ 3;  % m; length of problem domain
delta_t = 10^ -2;   % yr; time step (needs to be small, about 10^ -2 yr)
t_max = 10^ 1;      % yr; integration time of model 
M_0 = 4;            % m; net accumulation at upper end of problem domain
S_M = 4/ 10^ 4;     % 1/ m; net accumulation gradient (note change from 
                    % value given on wiki page, above)
bed_slope = 0.1;    % slope of glacier bed
Glen_n = 3;         % d'less; exponent in Glen's flow law
Glen_A = 10^ -16;   % Pa^ -3/ yr; coefficient in Glen's flow law
 
% Calculate derived constants.  
C = (2/ (Glen_n+ 2))* Glen_A* (rho* g)^ Glen_n; 
 
% Set up initial condition and plotting variables.  Initialize vectors.  
x = 0: delta_x: x_max;      % m; distances from upper end of domain for 
                            % each node
t = 0: delta_t: t_max;      % yr; times corresponding to each time step
z_bed = wrev(bed_slope* x); % m; bed height as a fn of x
num_x = numel(x);           % number of nodes
num_t = numel(t);           % number of time steps
D = zeros(1, num_x- 1);     % m^ 2/ yr; diffusivities at staggered points
H = zeros(1, num_x);        % m; thickness of ice at different points
M = zeros(1, num_x);        % m/ time step; net mass balance per time step
num_D = numel(D);           % number of points in offset grid
phi = zeros(1, num_D);      % m^ 2/ yr?; flux of matter between nodes
 
% Calculate the net mass balance per time step within the problem domain.  
M = (M_0- S_M* x)* delta_t; % m/ time step
 
% Step through time.  
for count1 = 1: 1: num_t; 
    for count2 = 1: 1: num_D; 
        D(count2) = C* ((H(count2)+ H(count2+ 1))/ 2)^ (Glen_n+ 2)* ...
            ((H(count2+ 1)+ z_bed(count2+ 1))- ...
            (H(count2)+ z_bed(count2)))^ (Glen_n- 1); 
        phi(count2) = D(count2)* (((H(count2+ 1)+ z_bed(count2+ 1))- ...
            (H(count2)+ z_bed(count2)))/ delta_x); 
    end
    for count3 = 2: 1: (num_x- 1); 
        H(count3) = H(count3)+ (phi(count3)- phi(count3- 1))/ ...
            delta_x+ M(count3); 
%         H(count3) = H(count3)- (phi(count3)- phi(count3- 1))/ ...
%             delta_x+ M(count3); 
    end
    H(1) = 0; 
    for count4 = 1: 1: num_x; 
        if H(count4) < 0; 
            H(count4) = 0; 
        end
    end
end
 
% Plot results.  
ice_top = z_bed+ H; % m; elevation of ice surface above lowest node
plot(x, ice_top, 'b-')
hold on
plot(x, z_bed, 'k-')
xlabel('Horizontal distance (m)')
ylabel('Elevation (m)')