f = [4 2; 
     4 3; 
     4 3; 
     4 4];
r = [1 0; 
     1 0; 
     1 -1; 
     0 0];

gamma = 0.9;

Q = rand(4,2);
%Q = 0.2*ones(4,2);
Q(end,:) = 0;
[m,n] = size(Q);

maxiter = 100;
epsilon = 0.001;

Qold = Q;
for i = 1 : maxiter,
    i, Q(1:3,:), pause;
    Q(:) = gamma*max(Q(f,:),[],2); Q = Q + r;
    if max(max(abs(Q-Qold))) < epsilon, break; end;
    Qold = Q;
end;

Q

return

for i = 1 : maxiter,
    i, Q(1:3,:), pause(.01);
    for x = 1 : 4,
        for u = 1 : 2,
            x1 = f(x,u);
            Q(x,u) = r(x,u) + gamma*max(Q(x1,:));
        end;
    end;
end;