I've been toying with procedural generation a lot recently but this is the first time I've come up with my own algorithm from scratch.

I'm very pleased with the results.

[edit: code updated]

set window on
maximize window
set display mode 1024,768,32
sync on
sync rate 0
randomize timer()
dim mapw(0)
dim maph(0)
w = 300
h = 300
mapw(0) = w
maph(0) = h
dim map(w,h)
maxres = 3^5
dim tmp(maxres,maxres)
create bitmap 1,w,h
set current bitmap 0
dim colour(512)
rem water
for i = 0 to 39
colour(i) = gradient(0,80,80, 20,140,140, i/39.0)
next i
rem beach
for i = 0 to 69
colour(i+40) = gradient(20,140,140, 255,255,140, (i^1.5)/(69^1.5))
next i
rem grass
for i = 0 to 109
colour(i+110) = gradient(255,255,140, 0,80,10, (i^.2)/(109^.2))
` colour(i+120) = rgb(110-i*.8,150-i*.7,30+i^.6)
next i
rem hills
for i = 0 to 74
colour(i+220) = rgb(40+i*1.6,40+i,i^1.1)
next i
rem snowy peaks
for i = 0 to 5
h = 55+i*40
colour(i+295) = rgb(h,h,h)
next i
colour(301) = 65535
iX = 5
iY = 5
size = 3^5 -1 :` size must be a power of 3 -1 or squares wont display properly.
res = 5
for x = 0 to 300
ink colour(x),0
box x*3,400,x*3+2,500
next x
do
rem clear map
boxdata(0,0,mapw(0),maph(0),0)
rem First step is to plonk a solid box on the map.
rem We will sculpt this into our island.
boxdata(iX,iY, iX+size-3, iY+size-3, 128)
rem Sculpt island.
island(iX,iY,size,res)
set current bitmap 1
for y = 1 to maph(0)-1
for x = 1 to mapw(0)-1
ink colour(map(x,y)),0
dot x,y
next x
next y
set current bitmap 0
copy bitmap 1,1,1,w-1,h-1, 0,41,41,40+w-1,40+h-1
rem monochrome heightmap.
set current bitmap 1
for y = 1 to maph(0)-1
for x = 1 to mapw(0)-1
c = int(map(x,y) * .85)*256
if c=0 then c=rgb(255,0,255)
ink c,0
dot x,y
next x
next y
copy bitmap 1,1,1,w-1,h-1, 0,441,41,440+w-1,40+h-1
sync:wait key
loop
end
function island(x,y,size,res)
if res > 1 then island(x,y,size,res-1)
scale = (size+res)/(3.0^res)
for v = 0 to 3^res-1
for u = 0 to 3^res-1
rem Get co-ordinates relative to scale.
x1 = x + u*scale
x0 = x1 - scale
x2 = x1 + scale
y1 = y + v*scale
y0 = y1 - scale
y2 = y1 + scale
rem Leave a border around the edge of the map.
if x0 <1 then x0 = 0
if y0 <1 then y0 = 0
if x2 >mapw(0) then x2 = mapw(0)
if y2 >maph(0) then y2 = maph(0)
rem Add up surrounding squares and apply weight to closest.
corners = map(x0,y0) + map(x2,y0) + map(x0,y2) + map(x2,y2)
sides = map(x1,y0) + map(x0,y1) + map(x2,y1) + map(x1,y2)
centre = map(x1,y1)
total = (corners*8 + sides*6 + centre*24)/res
rem random chance of land based on neighbours
tmp(u,v) = (total*res/60.0) - rnd(total/60.0)
if tmp(u,v) < 0 then tmp(u,v) = 0
if tmp(u,v) > 300 then tmp(u,v) = 300
next v
next u
rem write data to map
for v = 0 to 3^res-1
for u = 0 to 3^res-1
x0 = x+u*scale
y0 = y+v*scale
x1 = x+(u+1)*scale -1
y1 = y+(v+1)*scale -1
boxdata(x0,y0,x1,y1,tmp(u,v))
next u
next v
endfunction
`//
rem Fills a square block of data with the same value.
function boxdata(l,t,r,b,dat)
for y = t to b
for x = l to r
map(x,y) = dat
next x
next y
endfunction
`//
function gradient(r1,g1,b1,r2,g2,b2,pos#)
r = r1*(1-pos#) + r2*pos#
g = g1*(1-pos#) + g2*pos#
b = b1*(1-pos#) + b2*pos#
colour = rgb(r,g,b)
endfunction colour
`//