function R=ice1A(N,M,T,ic);
% ICE1A Solves simplified n=1 ice flow equation
%       h_t = 2/3 ( h_x (h-z0)^3 )_x
%    using a semi--implicit finite 
%    differences method.  Reference Mahaffy (1976).
%    dx=1/N; dt=T/M; R=2dt/(3 dx^2);
%    allows initial condition specified in ic;
%    Assume slab-on-a-slope:  z0 = 1-x.
%    Boundary conditions: h(t,0)=1, h(t,1)=0.
%    Displays 2 views: 3D and profile. Try:
%  "ice1A(10,50,0.1,[0 .5 1.0 .8 .4 0 0 0 0 0 0]+(1-(0:.1:1)))"
%  "ic = zeros(1,51)+(1-(0:.02:1));...
%    ic(20:30)=ic(20:30)+.2; ice1A(50,200,1.0,ic)"
% Ed Bueler 4/21/2000

dx=1/N; dt=T/M; x=0:dx:1; t=0:dt:T;

%allocate space for solution
w=zeros(M+1,N+1);
% insert initial condition
if not(all(size(ic) == [1 N+1]))
   R='Initial condition is of wrong size.'; return; end;
maxic=max(abs(ic))+1;
w(1,:)=ic;
w(1,1)=1; w(1,N+1)=0; %enforce boundary conditions at start

R=2*dt/(3* dx^2); Q=dt/dx; 
z=1-x(2:N);

for l=2:M+1
   w(l,1)=1; % boundary conditions
   w(l,N+1)=0; 
   % create tridiagonal N-1 X N-1 matrix A for left side
   wz=w(l-1,2:N)-z;
   aa = R * ( wz.^3 )'; % column vector
   A=spdiags([-aa ones(N-1,1)+2*aa -aa], -1:1, N-1, N-1);
   % b contains right side
   for k=2:N
      dw=w(l-1,k+1)-w(l-1,k-1);
      b(k-1)=w(l-1,k) + Q*wz(k-1)^2*dw*(dw/(2*dx)+1);
   end;
   b(1)=b(1)+aa(1)*w(l-1,1); % correct left boundary
   b(N-1)=b(N-1)+aa(N-1)*w(l-1,N+1); % correct right bdry
   w(l,2:N)=(A\b')';
   % compare: w(l,2:N)=tridiag(a,d,a,b,N-1);
   if any(abs(w(l,:))>5.0*maxic) break; end;
   % stops if size of solution blows up
end;

%display it
subplot(2,1,2), mesh(x,t,w), view(37.5,30);
xlabel('x axis'); ylabel('t axis'); zlabel('h axis');
title('3D view.  Slope defined by z0 shown.                    ');
% for profile plot, remove z0
for l=1:M+1
   w(l,:)=w(l,:)-[1 z 0];
end;
subplot(2,1,1), mesh(x,t,w), view(0,0);
xlabel('x axis'); zlabel('h axis'); 
title('n=1 slab-on-a-slope: slab thickness in profile.');