function T = func(U)         
N=size(U,1);
m=(N+1)/2;
H=1.0;         %Roughness factor. Vary between 0.5(very rough) and 1.5(very smooth)
n=log2(N-1);   %Number of iterations. Just solving N=(2^n)+1 for n.
p=2^(H*n);     %Scaling factor
T=U;
%Determine height in the middle of the square
T(m,m)=(U(1,1)+U(1,N)+U(N,1)+U(N,N))/4+(randn*p);  
%Determines height at midpoint of each edge of square
T(1,m)=(U(1,1)+U(1,N))/2+(randn*p);                    
T(m,N)=(U(1,N)+U(N,N))/2+(randn*p);                    
T(N,m)=(U(N,1)+U(N,N))/2+(randn*p);                    
T(m,1)=(U(1,1)+U(N,1))/2+(randn*p); 
%Reset the corners and run the function for each of the sub-squares created
if N > 3
    T(1:m,1:m)=func(T(1:m,1:m));
    T(1:m,m:N)=func(T(1:m,m:N));
    T(m:N,1:m)=func(T(m:N,1:m));
    T(m:N,m:N)=func(T(m:N,m:N));
end
end