function y = rbf_eval_stab(P,Fval,Pnt,rbf_func,varargin)
% Valutazione di una rbf con il trattamento dell'offset

% Dimensione del problema pari alle colonne di P
n = size(P,1);
dim = size(P,2);

% Inizializzazione della matrice K
K = zeros(n);

% Calcolo dell'offset
offset = [ones(n,1),P]\Fval;
if dim==1
    % Problema monodimensionale
    offsetP = offset(1) + offset(2)*P;
end
if dim==2
    % Problema monodimensionale
    offsetP = offset(1) + offset(2)*P(:,1) + offset(3)*P(:,2);
end


% Nuovi punti senza offset
Fval = Fval - offsetP;

% Ciclo sulle righe
for i = 1:n,
    % Ciclo sulle colonne
    for j = 1:n
        % Calcolo il valore di r_ij
        r = norm(P(i,:)-P(j,:));
        
        % Assegno il valore
        rbf_val = feval(rbf_func,r,varargin{:});
        K(i,j) = rbf_val;
    end
end

% Termine F
F = Fval;

% Risoluzione del sistema lineare
acoeff = K\F;


% Valutazione dell'interpolante
% -----------------------------


% Inizializzazione di y alla dimensione di Pnt
np = size(Pnt,1);
y = zeros(np,1);

% Ciclo su tutti i punti da valutare
for k=1:np
    % valutazoine di tutti i punti della rbf
    fval = zeros(length(acoeff),1);
    for i=1:length(acoeff)
        % valuto distanze rispetto al punto in considerazione
        r = norm(P(i,:)-Pnt(k,:));
        % valutazione della rbf
        rbf_val = feval(rbf_func,r,varargin{:});
        fval(i) = rbf_val;
    end
    
    % valutazione dell'interpolante nel punto
    y(k) = fval'*acoeff;
end

% aggiunta dell'offset
if dim==1
    % Problema monodimensionale
    offsetVal = offset(1) + offset(2)*Pnt;
end
if dim==2
    % Problema monodimensionale
    offsetVal = offset(1) + offset(2)*Pnt(:,1) + offset(3)*Pnt(:,2);
end

y = y + offsetVal;
