function Tnew = solveQuadratic(T, F, i, j, m, n)

% determine which derivative approximation is maximal - forwards or backwards
% i direction
if i == 1
    T1 = T(i+1, j, 1);
elseif i == m
    T1 = T(i-1, j, 1);
else
    T1 = min(T(i+1, j, 1), T(i-1, j, 1));
end
% j direction
if j == 1
    T2 = T(i, j+1, 1);
elseif j == n
    T2 = T(i, j-1, 1);
else
    T2 = min(T(i, j+1, 1), T(i, j-1, 1));
end

% solve
Tdiff = T1 - T2;
Fij = F(i, j);
if abs(Tdiff) < Fij
    % both i and j directions approximate the gradient - solve 2D quadratic equation
    Tnew = 0.5*(T1 + T2 + sqrt(2*Fij*Fij - Tdiff*Tdiff));
else
    % only i or j direction approximates the gradient - just add appropriate one to F
    Tnew = min(T1, T2) + Fij;
end 

