gamma = 0.25;
epsilon = 1e-12; 
t = 0.9;
% initialization of Q
for x = -1:1,
    for v = -1:1,
        Q(x+2,v+2) = f(x,v);
    end
end
% value iterations
k = 0;
while norm(Q-TQ(Q,gamma)) > epsilon,
    for x = -1:1,
        for v = -1:1,
            Q(x+2,v+2) = Q(x+2,v+2) - t*(Q(x+2,v+2) - ...
                TQlearning(Q,x,v,gamma));
            k = k+1;
        end
    end
end
Q
k

function Y = TQlearning(Q,x,v,gamma)
%
% Y = TQlearning(Q,x,v,gamma)
%
% evaluates the Bellman Q-operator at only (x,v)
%
for u = -1:1,
    temp(u+2) = gamma*Q(F(x,v)+2,u+2);
end
Y = f(x,v) + min(temp);
end

function Y = TQ(Q,gamma)
%
% Y = TQ(Q,gamma)
%
% evaluates the Bellman Q-operator
%
for x = -1:1,
    for v = -1:1,
        for u = -1:1,
            temp(u+2) = gamma*Q(F(x,v)+2,u+2);
        end
        Y(x+2,v+2) = f(x,v) + min(temp);
    end
end
end

function cost = f(x,u)
%
% cost = f(x,u)
%
% computes the incremental cost
%
cost = x^2 + u^2;
end

function xnew = F(x,u)
%
% xnew = F(x,u)
%
% computes the value of the new state given the current state x and the
% input u
%
if x == -1,
    if u == -1 | u == 0, 
        xnew = -1;
    else
        xnew = 0;
    end
elseif x == 0,
    if u == -1, 
        xnew = -1;
    elseif u == 0,
        xnew = 0;
    else
        xnew = 1;
    end
else
    if u == -1,
        xnew = 0;
    else 
        xnew = 1;
    end
end
end