I have experiemented with using matlabs built-in ode solvers to integrate a Shallow ice . Surprisingly the ode solvers are able to take much greater timesteps than what theoretical stability criteria would suggest. I am really impressed by how well it works although I will not recommend using this approach:
The main advantage of the ode solver approach is that the code can be made extremely short and readable. It could therefore serve as a nice introduction to the Shallow Ice Approximation.
Mail me if you would like a piece of the code. It is short and easy to read.
brown is the bedrock. It is a modification of matlabs peaks function.
transparent white is the glacier surface.