clear; 
load B 
load L 
load F 
S = size(F); 
dim = S(1); 
R = zeros(dim,dim); 
epsilon  = .0001; 
maxiterate = 100; 



i = 1;  
D(1) = 1000; 
while(D(i) > epsilon) 
    i = i + 1; 
    old = R; 
    R = -R * R * B * (L)^(-1) - F * (L)^(-1); 
    diff = 0; 
    for k = 1:dim 
        for l = 1:dim 
            diff = max(diff, abs(old(k,l) - R(k,l))); 
        end 
    end 
    D(i) = diff; 
end 

save R R