## The first set of code plays only one set of games for the
## gambler's ruin problem. It is short and thus easy to de-bug.
## After we have one set of code working, it is easier to modify
## it to do a long loop of sets of games. 


p:= 0.5: # initializes the probability of heads = probability you gain $1
x := 10: # initializes the amount of money you have
y := 10: # intitializes the amount of money your opponent has
t := 0:  # intializes the number of coin flips that have occurred
broke := 0: # initializes whether you are broke or not 
## The following code plays one set of games - 
  ##  until you are broke or your opponent is broke ## 
while x>0 and x<x+y do   
u := rand()/10^12;  ##pick a random number in (0,1)
if (u<p) then       ## determine the outcome of the coin toss
x := x + 1;         ## if it's heads, you gain and 
y := y - 1;         ## your opponent loses one dollar
t := t + 1;         ## and you have flipped another coin
else                ## else it is the reverse
x := x - 1;
y := y + 1;
t := t + 1;
end if;            ## remember to end your ifs
if (x = 0) then broke := 1 end if;   ## record if you are broke
end do:            ## remember to end your dos


#####################################################################


## The following code modifies the above code in order to run a 
## single long simulation. 


number_simulations := 100:
t := Array(1..number_simulations):  ## array to store number of flips
broke := 0: # initializes whether you are broke or not 
## The following code plays one set of games - 
  ##  until you are broke or your opponent is broke ## 
for i from 1 to number_simulations by 1 do
p:= 0.5: # initializes the probability of heads = probability you gain $1
x := 10: # initializes the amount of money you have
y := 10: # intitializes the amount of money your opponent has
while x>0 and x<x+y do   
u := rand()/10^12;  ##pick a random number in (0,1)
if (u<p) then       ## determine the outcome of the coin toss
x := x + 1;         ## if it's heads, you gain and 
y := y - 1;         ## your opponent loses one dollar
t[i] := t[i] + 1;         ## and you have flipped another coin
else                ## else it is the reverse
x := x - 1;
y := y + 1;
t[i] := t[i] + 1;
end if;            ## remember to end your ifs
## record number of times you are broke
if (x = 0) then broke := broke + 1 end if;   
end do:            ## remember to end your dos
end do:            ## remember to end ALL of your dos



#####################################################################

## The following code modifies the above code to vary x for y = 10



number_simulations := 100:
t := Array(1..20, 1..number_simulations):  ## array to store number of flips
broke := Array(1..20): # initializes whether you are broke or not 
## The following code plays one set of games - 
  ##  until you are broke or your opponent is broke ## 
for j from 1 to 20 by 1 do
for i from 1 to number_simulations by 1 do
p:= 0.5: # initializes the probability of heads = probability you gain $1
x := j:  # intializes the amount of money you have
y := 10: # intitializes the amount of money your opponent has
while x>0 and x<x+y do   
u := rand()/10^12;  ##pick a random number in (0,1)
if (u<p) then       ## determine the outcome of the coin toss
x := x + 1;         ## if it's heads, you gain and 
y := y - 1;         ## your opponent loses one dollar
t[j,i] := t[j,i] + 1;         ## and you have flipped another coin
else                ## else it is the reverse
x := x - 1;
y := y + 1;
t[j,i] := t[j,i] + 1;
end if;            ## remember to end your ifs
## record number of times you are broke
if (x = 0) then broke[j] := broke[j] + 1 end if;   
end do:            ## remember to end your dos
end do:            ## remember to end ALL of your dos
end do:            ## remember to end ALL of your dos



## The following code finds the averages from above and plots them against 
## the amount of money you started with. Remember y = 10 and p = .5



average_time := Array(1..20):
for j from 1 to 20 by 1 do
for i from 1 to number_simulations by 1 do
average_time[j] := average_time[j] + t[j,i];
end do:
end do:
average_time := evalf(average_time/number_simulations);
broke := evalf(broke/number_simulations);


with(plots):
points:= {seq([i,average_time[i]],i=1..20)}:
pointplot(points);