mmax = 100;

for m=1:mmax,
 for n=1:m,
 l2 = m*m + n*n;
 l = round(sqrt(l2));
 if (l*l-l2 == 0)
  fprintf (1,'%d*%d = %d*%d + %d*%d \n',l,l,m,m,n,n);
  for mp = 1:sqrt(l)-1,
  np = round(sqrt(l-mp*mp));
  if (np*np + mp*mp == l) 
   fprintf (1,'%d = %d*%d + %d*%d \n',l,mp,mp,np,np);
  end;
  end;
 end;
end; 
end;
