function S = Sfun2D(b)
% computation of an error function for an ODE model
% INPUT: b - vector of parameters

global ndata xdata x0

%% discrete model 
% (nested function, uses parameters b(1) and b(2) of the main function)
% input model as Dx_n = f(x_n) where Dx_n = x_n+1 - x_n)
    function dx = f(x)
        dx = zeros(2,1);
        dx(1) = -b(1)*x(1)*x(2);
        dx(2) = b(1)*x(1)*x(2) - b(2)*x(2);
    end

%% model integration

nsol = [0:max(ndata)];
x = zeros(length(nsol)+1,length(x0));
% note the index change: x(0) represents the value x at n = 0;
xsol(1,:) = x0;  
for n = nsol(1:end-1)+1
    xsol(n+1,:) = f(xsol(n,:))' + xsol(n,:) ;
end
    

%% plot result of the integration

figure(1)
for i = 1:2
    subplot(1,2,i)
    plot(ndata,xdata(:,i),'x','MarkerSize',10);
    hold on
    plot(nsol,xsol(:,i),'o');
    hold off
    ylabel(['x(' num2str(i) ')']);
end
drawnow

%% find predicted values x(tdata)

xpred = interp1(nsol,xsol,ndata);

%% compute total error

S = 0;
for i = 1:length(ndata)
    S = S + sum((xpred(i,:)-xdata(i,:)).^2);
end

end