function[R] = squareradius(maxrad) 
%MATLAB program that calculates the first maxrad terms of this sequence: 
%R(r) is the minimal number of additions needed so that site (0,r) topples 
%in the abelian sandpile growth model with h = 2. 
M = 2.0*ones(3,3); 
% M is the number of grains per site. 
L = [0,1,0;1,-4,1;0,1,0]; 
% effect of a toppling 
id = [0,0,0;0,1,0;0,0,0]; 
%no effect 
R = [2]; 
%the first term 
add = 0; r = 1; breedte = 3; while r<maxrad M(r+1,r+1) = M(r+1,r+1)+1; 
%make an addition to the center of M 
add = add+1; while max(max(M))>3 
% stabilize M by topplings 
instab = 1.0*(M>3); if ~isequal(instab(:,1),zeros(breedte,1)) 
%r increases 
dM = imfilter(instab,L,'full'); 
%function from image processing toolbox 
M = imfilter(M,id,2,'full'); 
breedte = breedte + 2; r = r+1; 
R(r) = add; else dM = imfilter(instab,L); end M = M + dM; end end;