#sparse matrices - 3 diag matrix
1==1;

function L=lap1s(n)
#creates sparse discrete minus 1d laplacian on
#uniform grid scaled by h^2 i.e. 3 diag matrix with 
#2 on the main diag, -1 on 1st super and sub diagonals
#n size of L nxn\
#OUTPUT 
#L 3diag matrix
#
#Example: L=lap1s(4); full(L) 
 S=-ones(1,3*n-2);
 S(1:n)*=-2;
 I=[1:n,1:n-1,2:n];#row indices: diagonal; superdiag and subdiag
 J=[1:n,2:n,1:n-1];#column indices: diagonal; superdiag and subdiag
 L=sparse(I,J,S,n,n);
endfunction


function [u,x]=mixsolver(f,s,a,b,N)
#solver for mixed problem -u''=f 
#u(a)=s u'(b)=0 by FDM on equididstant mesh
#u_0=s Dirichlet bnd
# -u_{k-1}+2u_k -u_{k+1}=h^2*f(x_k) k=1..N-1
#u_N-u_{N-1}=h*0 Neumann bnd
#INPUT
#f - rhs function of ODE - vector implemented
#s - Dirichlet bnd
#a,b, - interval ends
#N - no of subintervals - N+1 no of meshpoints 

x=linspace(a,b,N+1); x=x';
F=f(x);
A=sparse(N+1,N+1);
A=lap1s(N+1,N+1);
h=(b-a)/N;
#we modify A and F:
F*=h*h;
#Dirichlet bnd at the left end
A(1,1)=1;A(1,2)=0;
F(1)=s;
#Neumann bnd at the right end - zero 
A(N+1,N)=-1;
A(N+1,N+1)=1; 
F(N+1)=h*0; 
u=A\F;
endfunction


function y=objectf(h,u,x)
#computes approximation of 
#\int_{x_1}^{x_N} |h-u|^3 by trpezoidal rule
#INPUT
#x- meshpoints
#h,u - values of h nd u at the meshpoints
#OUTPUT
#y=approximation of th square of integral by trapezoidal rule
di=h-u;
y=trapz(x,di.*di);

end



function y=realobf(s)
#objective function with one parameter
global f h a b n
[u,x]=mixsolver(f,s,a,b,2*n);
ind=n+1:2*n+1;
xx=x(ind);
hh=h(xx);
y=objectf(hh,u(ind),xx);
endfunction


#testing for f=0;s=1; h=0.5 then u=1 and \int_1^2 |u-h|^2=0.25
global f h a b n
f=@(x) 0*x;
h=@(x) 0.5 + 0*x;
a=0;b=2;
s=1;n=20;
y=realobf(s)

#optimization
s0=fminbnd(@realobf, -1,1)

#is s0==0.5?
er=s0-0.5 #OK!

f=@(x) (x>0.5).*(x<1.5);
h=@(x) 1-(2-x).*(2-x);
ss=linspace(0,1.5,50);
for k=1:length(ss),
 fs(k)=realobf(ss(k));
endfor
plot(ss,fs);pause(1);
s0=fminbnd(@realobf, -1,3)
[u,x]=mixsolver(f,s0,a,b,2*n);
[u1,x]=mixsolver(f,s0+0.1,a,b,2*n);
plot(x,u,";s0;",x,u1,";s=s0+0.1;",x,h(x),";h;")

