function simulation(Tmax, Xmax)
    tic
    % Run turtle pattern simulation
    % if Xmax is not set we use a default of 40
    % if even Tmax is not set we use a default of 1800

    if nargin < 2, Xmax = 40; end;   % Default max x value
    if nargin < 1, Tmax = 1800; end; % Defalt max t value

    % Set up constants for periodic pattern:
    dt = 0.1;			dx = 1;     % Step-sizes
    au = 0.02;	av = 0.03;          % Decay constants
    bu = 0.001;	bv = 0;             % Base expression rates
    Du = 0.01;	Dv = 0.15;           % Diffusion constants

    % Define the domain:
    Nt = floor(Tmax/dt) + 1;
    Nx = floor(Xmax/dx) + 1;
    T = dt*(0:(Nt-1));
    X = dx*(0:(Nx-1));

    % Set up slightly noisy initial conditions:
    noise = 0.005;
    lu = au*(1 + noise*(2*rand(Nx +2,Nx +2)-1));
    u = ones(Nx,Nx);
    v = u;

    % Reaction rules:
    f = @(u,v) lu.*(u.^2 + bu)./(v) - au*u;
    g = @(u,v) lu.*u.^2 - av*v + bv;

    % Display initial state:
    figure(1); hax1 = gca();    % create window 1 and draw u
    draw(hax1, u);
    figure(2); hax2 = gca();    % create window 2 and draw v
    draw(hax2, v);

    % Simulation loop:
    for t = 1:Nt -1

        % Expand u and v to use zero-flux boundary conditions
        % Note: this just works as our laplace2d applies the matrices
        % in the correct order
        ue = zeroFlux(u);
        ve = zeroFlux(v);

        % Calculate interdots (rk)s
        ue2 = ue + dt/2 * (f(ue, ve) + Du*laplace2d(ue, dx));
        ve2 = ve + dt/2 * (g(ue, ve) + Dv*laplace2d(ve, dx));

        % Calculate new u and v using the interdots
        ue = ue + dt * (f(ue2, ve2) + Du*laplace2d(ue2, dx));
        ve = ve + dt * (g(ue2, ve2) + Dv*laplace2d(ve2, dx));

        % Shrink u and v back to real size
        u = ue(1:Nx,1:Nx);
        v = ve(1:Nx,1:Nx);

        % Fixed Boundary Condition
%         uBound = 0;
%         vBound = 1; % vBound must not be 0
%         u = fixedBoundary(u, uBound);
%         v = fixedBoundary(v, vBound);

        % Display:
        if mod(t,50) == 0
	    draw(hax1,u);   % This is necessary if we want animation (calls drawnow)
	    draw(hax2,v);   % This is necessary if we want animation (calls drawnow)
       end;
    end;
    toc
end