function [e,J]=nnfun(x,Phi,Y,Nh)
% Objective function and Jacobian for neural network example

N=size(Y,1);     % Number of data points
Ni=size(Phi,2);  % Number of input neurons

% Extract weights and biases from x.
Wh=zeros(Ni,Nh);
Wh(:)=x(1:Ni*Nh);
Bh=zeros(1,Nh);
Bh(:)=x(Ni*Nh+1:Ni*Nh+Nh);
Wo=zeros(Nh,1);
Wo(:)=x(Ni*Nh+Nh+1:Ni*Nh+2*Nh);
Bo=x(Ni*Nh+2*Nh+1);

% Compute the output of the neural network.
V=Phi*Wh+ones(N,1)*Bh;
Ye=tanh(V)*Wo+ones(N,1)*Bo;

e=(Y-Ye);
plot([Y Ye e]); % Display the training.
drawnow;

% Compute the Jacobian if necessary.
if ( nargout > 1 )
   thV=tanh(V');
   dBo=ones(1,N);
   dWo=thV;
   dBh=(Wo*ones(1,N)).*(ones(Nh,N)-(thV).^2);
   dWh=zeros(Ni*Nh,N);
   for i=1:Nh,
      dWh(Ni*(i-1)+1:Ni*i,:)=(ones(Ni,1)*dBh(i,:)).*Phi';
   end;
   J=-[dWh; dBh; dWo; dBo]';
end;
