Purpose of writing this program:
I started reading RichardDawkins.net, read a few of his books, and later started reading Pharyngula. "Climbing Mount Improbable" and "The Selfish Gene" are a better description of evolution than anything else I've ever seen. I'll admit to having some doubts as to how an unintelligent evolutionary process could produce complexity, but they were demolished by those books. At some point one of those sites linked to something that linked to some creationist site (I can't find it now) claiming that a computer program reproducing eye evolution referenced in "Climbing Mount Improbable" doesn't actually exist. I found that paper, and it in fact describes the small changes that could occur producing advanced eyes from simple light detecting cells in order to calculate the time necessary, and does not claim the authors produced a computer program reproducing this.
I wrote a quick and dirty version of such a program.

Why I am posting this:
I wanted this code to be openly available. This is intended mostly for engineering students, or people in a university setting in general. Apparently there are more creationists among engineers and computer scientists (and physicians and mathematicians) than biologists, chemists, or physicists. I've met some of them. This program is written in MATLAB, which is an industry standard for engineering, and should be available to every engineering student in the country (probably the world). MATLAB is a relatively easy language, and the code is heavily commented, so I expect any engineer to be able to read it. Those of you who know creationist engineering students can challenge them to look at and run this program, and find a problem.

Preemption of criticism:
I know there are issues with the ray tracer: I made assumptions when writing it that the algorithm naturally did not respect (damn unguided evolution). The most obvious is that it is two dimensional, and has constructed things that would never occur in three dimensions. Similarly, I am not a biologist, so nothing here represents the real development of an eye, nor is it subject to the same physical limitations. Finally, the fitness functions used are constant, and naturally don't take into account everything that would matter in an organism.
The program does demonstrate the evolution of complex structures that perform useful functions from simpler ones, which was MY goal in the first place. The inaccuracies mentioned don't change the fact that an unintelligent algorithm (anyone feel like arguing that?) creates complexity (same question) based on Darwinian principles of variation generated by random mutation selected by relative fitness. It just happens to be part of an artificial structure selected in an artificial universe according to its artificial laws of physics.
I know there are better ways to implement many parts of this. Maybe someday I will. If you modify something here and get a useful result, please send it to me, so I can post it here.
My next post will show how to use the program and some results and conclusions.

=======================================================================================================
%version 5
function [besteye,bestem,score,out]=evolveeye(n,emin,s)%data fields on emin: refractive index, x delta, y delta for each cell
%allows an eye to evolve from a flat plate of light sensing cells
format short g
if exist('s','var')
    rand('state',s);
end
embryo{1}=emin; ss=1; fn=1;                 %Set save to file parameters.
for tt=1:n                                  %generational loop:
    layer{1}=groweye(embryo{1});            %Create
    [score,out]=visiontest(layer{1});       %and test parent eye.
    best=[1,score];                         %Initialize structure to track most successful mutant.
    for ii=2:10                             %population loop:
        embryo{ii}=mutateall(embryo{1});    %duh
        layer{ii}=groweye(embryo{ii});      %Create
        [score,out]=visiontest(layer{ii});  %and test mutant eye.
        if score<best(2)                    %New best variant found,
            best=[ii,score];                %save appropriate data.
        end
    end
    [tt,best(2)]                            %Display progress.
    embryo{1}=embryo{best(1)};              %Best variant from this generation
    layer{1}=layer{best(1)};                %becomes parent to the next.
    [nm,nn,nl]=size(layer{1});
    c=reshape(layer{1}(:,1,:),nm,nl);       %Prepare shading for plot.
    x=reshape(layer{1}(:,2,:),nm,nl);       %Prepare x coords for plot.
    y=reshape(layer{1}(:,3,:),nm,nl);       %Prepare y coords for plot.
    figure(1);hold off;colormap('gray');contourf(x,y,c);%Shade refractive index.
    hold on;plot(x,y,'.-');%axis equal      %Plot current eye.
    if mod(tt,10)==0||tt==n                 %Store every tenth generation
        saveeye{ss,1}=layer{1};
        saveeye{ss,2}=embryo{1};
        ss=ss+1;
        if ss==11||tt==n                    %and save to file every hundreth generation.
            eval(sprintf('save eye%03d saveeye',fn));
            ss=1; fn=fn+1;
        end
    end
    pause(.01);                             %Necessary to update display.
end
score=best(2);                              %Set useful outputs.
besteye=layer{1}; bestem=embryo{1};

function layer=groweye(em)
%grows an eye from a "stem cell" defined by a location and property
%using programmed embryology defining change in location and
%properties for each new cell.
[nm,nn,nl]=size(em); layer=zeros(2*nm-1,nn,nl);
for ll=1:nl                                     %loop over all layers of eye:
    for mm=1:nm                                 %loop for radial position of eye:
        if ll==1&&mm==1                         %Start point is
            layer(nm,:,1)=em(1,:,1);            %defined by stem cell.
        elseif mm==1                            %First axial position is
            layer(nm,:,ll)=layer(nm,:,ll-1)+... %modified from cell below
                em(1,:,ll);                     %by embrylogy instructions.
        else                                    %Other cells are
            layer(nm+mm-1,:,ll)=...             %modified from next cell axially
                layer(nm+mm-2,:,ll)+em(mm,:,ll);%by embrylogy instructions
            layer(nm-mm+1,:,ll)=...             %in both directions horizontally.
                layer(nm-mm+2,:,ll)+[1,-1,1].*em(mm,:,ll);
        end
    end                                         %end radial loop
end                                             %end layer loop

function em2=mutateall(em)
%Add random data to embryology and correct to physical posibilities.
[nm,nn,nl]=size(em); mutmag=.001; sizefreq=.0001;
%mutmag=.0001 goes nowhere after 8000 iterations
em2=em; resizer=rand; r2=rand;                  %Random value to add new axial or radial positions.
                                                %Matlab command provides
                                                %equal probability distribution between 0 and 1.

if resizer<sizefreq                             %sizefreq% chance of
    tmpi=floor(nm*r2)+2;                        %inserting new layer
    if tmpi-1>0&&tmpi<nl                        %at random position
        em2(:,:,tmpi+1:end+1)=em(:,:,tmpi:end); %either by interpolating
        em2(:,:,tmpi)=(em2(:,:,tmpi+1)+em2(:,:,tmpi-1))/2;
    else
        em2=em; em2(:,:,end+1)=em(:,:,end);     %or copying onto the outside.
    end
elseif nl>2&&resizer<2*sizefreq                 %given enough layers,
    tmpi=floor(nl*r2)+1;                        %sizefreq% chance of
    em2=em(:,:,[1:tmpi-1,tmpi+1:end]);          %losing random layer.
elseif resizer<3*sizefreq                       %sizefreq% chance of
    tmpi=floor((nm+1)*r2)+1;                    %randomly inserting new radial station
    if tmpi>1&&tmpi<nm                          %by either interpolation,
        em2(tmpi+1:end+1,:,:)=em(tmpi:end,:,:);
        em2(tmpi,:,:)=(em2(tmpi+1,:,:)+em2(tmpi-1,:,:))/2;
    elseif tmpi==1%add center
        em2(tmpi+1:end+1,:,:)=em(tmpi:end,:,:); %copying a new center,
        em2(1,:,:)=em(1,:,:);
    else%add end
        em2=em; em2(end+1,:,:)=em(end,:,:);     %or copying the outermost location.
    end
elseif nm>2&&resizer<4*sizefreq                 %given enough radial stations,
    tmpi=floor(nm*r2)+1;                        %sizefreq% chance of
    em2=em([1:tmpi-1,tmpi+1:end],:,:);          %losing random radial station.
end
[nm,nn,nl]=size(em2);
em2=em2+mutmag*randn(nm,nn,nl);                 %Generate variation in genome,
em2=fixeye(em2);                                %and repair any impossibilities.

function em2=mutate(em)
%Add random data to embryology and correct to physical posibilities.
[nm,nn,nl]=size(em); mutmag=.001; sizefreq=.0025;
em2=em; resizer=rand; r2=rand;                  %Random value to add new axial or radial positions.
                                                %Matlab command provides
                                                %equal probability distribution between 0 and 1.

if resizer<sizefreq                             %sizefreq% chance of
    tmpi=floor(nm*r2)+2;                        %inserting new layer
    if tmpi-1>0&&tmpi<nl                        %at random position
        em2(:,:,tmpi+1:end+1)=em(:,:,tmpi:end); %either by interpolating
        em2(:,:,tmpi)=(em2(:,:,tmpi+1)+em2(:,:,tmpi-1))/2;
    else
        em2=em; em2(:,:,end+1)=em(:,:,end);     %or copying onto the outside.
    end
elseif nl>2&&resizer<2*sizefreq                 %given enough layers,
    tmpi=floor(nl*r2)+1;                        %sizefreq% chance of
    em2=em(:,:,[1:tmpi-1,tmpi+1:end]);          %losing random layer.
elseif resizer<3*sizefreq                       %sizefreq% chance of
    tmpi=floor((nm+1)*r2)+1;                    %randomly inserting new radial station
    if tmpi>1&&tmpi<nm                          %by either interpolation,
        em2(tmpi+1:end+1,:,:)=em(tmpi:end,:,:);
        em2(tmpi,:,:)=(em2(tmpi+1,:,:)+em2(tmpi-1,:,:))/2;
    elseif tmpi==1%add center
        em2(tmpi+1:end+1,:,:)=em(tmpi:end,:,:); %copying a new center,
        em2(1,:,:)=em(1,:,:);
    else%add end
        em2=em; em2(end+1,:,:)=em(end,:,:);     %or copying the outermost location.
    end
elseif nm>2&&resizer<4*sizefreq                 %given enough radial stations,
    tmpi=floor(nm*r2)+1;                        %sizefreq% chance of
    em2=em([1:tmpi-1,tmpi+1:end],:,:);          %losing random radial station.
end
[nm,nn,nl]=size(em2);
c=[ceil(nm*rand),ceil(nn*rand),ceil(nl*rand)];  %At random location in genome,
em2(c(1),c(2),c(3))=em2(c(1),c(2),c(3))+mutmag*randn;%generate point mutation,
em2=fixeye(em2);                                %and repair any impossibilities.

function em2=fixeye(em)
%Add random data to embryology and correct to physical posibilities.
constr=.2;
em2=em;
[nm,nn,nl]=size(em2);
test=groweye(em2);                             %Make eye to look for problems.

% c=reshape(test(:,1,:),2*nm-1,nl);       %Make data and show running plot.
% x=reshape(test(:,2,:),2*nm-1,nl);       %Make data and show running plot.
% y=reshape(test(:,3,:),2*nm-1,nl);
% figure(1);hold off;colormap('gray');contourf(x,y,c);
% hold on;plot(x,y,'.-');%axis equal
km=nm-1;                                    %Embryo to eye index conversion
for ll=1:nl                                 %Loop over all layers of eye:
    for mm=1:nm                                %Loop for radial position of eye:
        if mm==1                               %Center of layer is constrained to axis of eye.
            em2(mm,2,ll)=0;
        elseif em2(mm,2,ll)<0                 %Radial distance between points
            em2(mm,2,ll)=0.01;                 %is constrained by minimum value.
        elseif em2(mm,2,ll)>constr
            em2(mm,2,ll)=constr;
        end
        if mm==1&&ll==1
            em2(mm,3,ll)=0;
        elseif mm==1&&em2(mm,3,ll)<0
            em2(mm,3,ll)=0;
        elseif abs(em2(mm,3,ll))>constr        %Radial distance between points
            em2(mm,3,ll)=sign(em2(mm,3,ll))*constr;             %is constrained by minimum value.
        end
        if test(mm+km,1,ll)<1               %If refractive index is not physically
            em2(mm,1,ll)=em2(mm,1,ll)+1-test(mm+km,1,ll);%possible, correct this.
        end
        test=groweye(em2);
        if ll>1                               %For layers above the retina,
            if mm==1||km<=1                   %y location of layer below is calculated
                y=test(mm+km,3,ll-1);          %knowing that the center of that layer
            else                               %is directly below this center,
                tmpi=find(test(nm:end,2,ll-1)>test(mm+km,2,ll));
                if length(tmpi)==0             %or by finding the points under the current point
                    tmpi=nm;
                end
                tmpi=tmpi(1)+km;              %and interpolating the height below the current point.
               

y=(test(mm+km,2,ll)-test(tmpi-1,2,ll-1))/(test(tmpi,2,ll-1)-test(tmpi-1,2,ll-1))*(test(tmpi,3,ll-1)-test(tmpi-1,3,ll-1))+test(

tmpi-1,3,ll-1);
            end
            if test(mm+km,3,ll)<y             %The current layer is constrained to not be below
                em2(mm,3,ll)=em2(mm,3,ll)+y+0.01-test(mm+km,3,ll);%the underlying layer.
            elseif test(mm+km,3,ll)>y+constr       %The current layer is constrained to not be >constr
                em2(mm,3,ll)=em2(mm,3,ll)+y+constr-test(mm+km,3,ll);%over the underlying layer.
            end
        end
        test=groweye(em2);                   %Test eye is regenerated to account for modifications.
    end                                     %end radial looptest=groweye(em2);
end
% c=reshape(test(:,1,:),2*nm-1,nl);       %Make data and show running plot.
% x=reshape(test(:,2,:),2*nm-1,nl);       %Make data and show running plot.
% y=reshape(test(:,3,:),2*nm-1,nl);
% figure(2);hold off;colormap('gray');contourf(x,y,c);
% hold on;plot(x,y,'.-');%axis equal

Continued here.