
maxIterations = 500;
gridSize = 1000;
rlim = [-0.748766713922161, -0.748766707771757];
ilim = [ 0.123640844894862,  0.123640851045266];


% Setup
t = tic();
x = gpuArray.linspace( rlim(1), rlim(2), gridSize );
y = gpuArray.linspace( ilim(1), ilim(2), gridSize );
[xGrid,yGrid] = meshgrid( x, y );
z0 = complex( xGrid, yGrid );
count = gpuArray.zeros( size(z0) );

% Calculate
z = z0;
for n = 0:maxIterations
    z = z.*z + z0;
    inside = abs( z )<=2;
    count = count + inside;
end
count = log( count+1 );

% Show
naiveGPUTime = toc( t );
%imagesc( x, y, count )
%axis image
%title( sprintf( '%1.2fsecs (naive GPU) = %1.1fx faster', ...
%    naiveGPUTime, cpuTime/naiveGPUTime ) )

fprintf( '%1.2fsecs (with GPU but naively)\n', naiveGPUTime );

