
disp("optimal contro; find s=arg min_s \int_0.5^1 |u-H|^2 dx with H prescribed");
disp("and u solving -u''=f; u(0)=s; u'(1)=0");
graphics_toolkit("gnuplot");

0==0;
%FDM for a linear BVP -u''+cu=f in (a,b) u(a)=alpha u'[(b)=beta
f=@(x) -2*ones(size(x));
N=10;s=1;
[u,x]=fdmbvpsolmix(f,0,1,2*N,0,s,0);

function [y,u,x]=Error(s=1,f=@(x) -2*ones(size(x)),H=@(x) 20*ones(size(x)),N=20)
   [u,x]=fdmbvpsolmix(f,0,1,2*N,0,s,0);
   umH=u-H(x);
   ers=umH.*umH;
   h=x(2)-x(1);% [x(N+1),x(2*N+1)]
   y=trapz(ers(N+1:2*N+1))*h;
end   

function y=F(s)
  global ff HH uu xx NN
  NN=20;
  [y,uu,xx]=Error(s,ff,HH,NN);
end  


  global ff HH uu xx NN
  ff=@(x) -2*ones(size(x));
  HH=@(x) 20*ones(size(x))
  
  ss=linspace(15,25);
  k=1;
  for s=ss,
    FF(k++)=F(s);
  end
  s0=fminbnd(@F,15,25)
  plot(ss,FF,s0,Fs0=F(s0));pause(2);  
  Fs0
  I=NN:2*NN+1;
  plot(xx(I),uu(I),xx(I),HH(xx(I)))