DECLARE EXTERNAL FUNCTION poisson
OPTION BASE 0
RANDOMIZE 
DIM a(100)

FOR k=1 TO 10000
   LET x=RND
   LET j=0
   LET s=poisson(0)
   DO WHILE s<x
      LET j=j+1
      LET s=s+poisson(j)
   LOOP
   LET a(j)=a(j)+1
NEXT k

SET WINDOW -3,31,-100,2000
DRAW axes(1,100)
FOR j=0 TO 30
   PRINT j,a(j)
   PLOT LINES:j-0.5,0;j-0.5,a(j);j+0.5,a(j);j+0.5,0;j-0.5,0
NEXT j

END

EXTERNAL FUNCTION poisson(k)
LET lambda=10
IF k=0 THEN
   LET poisson=EXP(-lambda)
ELSE
   LET t=1
   FOR j=1 TO k
      LET t=t*j
   NEXT j
   LET poisson=EXP(-lambda)*lambda^k/t
END if
END FUNCTION