$M=5000; # granularity $N=2000; # number of terms in the sum $lambda=1; # must be <= 1 $alpha=1; $beta=1/$alpha; open(OUT,">solve.txt"); for ($k=0; $k<=$M; $k++) { $a[$k]=1; } for ($iter=0; $iter<10; $iter++) { $sum=0; for ($k=0; $k<=$M; $k++) { $b[$k]=0; $x=$k/$M; for ($n=1; $n<=$N; $n++) { $y=($lambda/($x+$n))**$beta; $y=int(0.5+ $M*$y); ## must be: 0 <= y < M if ($y<=$M) { $b[$k]+=$a[$y]*$beta*($lambda**$beta)*(($x+$n)**(-1-$beta)); } else { print "ERROR\n"; # ..... error? exit? } } $sum+=$b[$k]; if ($k%20 ==0) {print OUT "$iter\t$x\t$b[$k]\n"; } } for ($k=0; $k<=$M; $k++) { $a[$k]=$b[$k]*$M/$sum; } print "$iter | $M | $a[0] | $a[$M] | $sum\n"; } close(OUT);