

%solving nonlinear problem defined by IVP i.e.
%for given IVP y'=f(x,y) y(x0)=y0 find y0 s.t. y(x1)=y1 (x1>x0)
%two solutions: formulate it as a nonlinear univariate problem and 
%solve it using fzero/fsolve - nonlinear function may use lsode()
%Second solution: integrate backward the problem from x1 t ox0 
%test for y'=y*(1-y)*sin(y*y+x) y(0)=s y(1)=0.5
global data
data.f=@(x,t) x*(1-x)*sin(x*x+t);
data.t0=0;
data.t1=1;
data.x1=0.5;
function y=G(s)
%function returning the solution at t1
%of IVP x'=f(t,x) x(t0)=s   
%f - functon handle t0,t1 parameters in structure which is global variable
global data
t1=data.t1;
t0=data.t0;
f=data.f;
X=lsode(f,s,[t0,t1]);
y=X(2,:)';
end

[G(0),G(1)]-[0,1]


function y=Gn(s)
global data
y=G(s)-data.x1;
end  

[s0,xt1]=fzero(@Gn,[0,1])
t1=data.t1
t0=data.t0
x1=data.x1
f=data.f;
t=linspace(0,1.5);
X=lsode(f,s0,t);
plot(t,X,t,x1*ones(size(t)));pause(1);pause(1);

%we integrate backward 
T=data.t1-data.t0;
%s=t0-t =-t dx/ds=dx/dt(t(s)(dt/ds)=-f(x,-s)
g=@(x,s) -f(x,t1-s)
X=lsode(g,0.5,[0,T]);s1=X(2);
X=lsode(f,s1,[0,1]);
X(2)-x1
s0-s1%==0?

%boundary value problem BVP
%how to solve u''+u*u=f (f=exp(t) or exp(-t) or t+1
%u(0)=0 u(1)=0
%shooting method let 
%u(t;s) be the solution of IVP for this ODE with u(0)=0 u'(0)=s
%so we have to solve nonlinear problem g(s)=0 with g:=u(1;s)

%defining g(s)
global f
f=@(t) exp(t);
function y=g(s)
  global f
  pw=@(x,t) [x(2);-x(1)*x(1)+f(t)];
  X=lsode(pw,[0;s],[0,1]);
  y=X(2,1);
end
[s0,gs0]=fsolve(@g,0)  
t=linspace(0,1);
X=lsode(@(x,t) [x(2);-x(1)*x(1)+f(t)],[0;s0],t);
%plot(t,x=X(:,1));
[min(x),max(x)]
[s1,gs1]=fsolve(@g,20)  
X=lsode(@(x,t) [x(2);-x(1)*x(1)+f(t)],[0;s1],t);
plot(t,x,[";s=",num2str(s0),";"],t,x1=X(:,1),[";s=",num2str(s1),";"]);
pause(1);
[min(x1),max(x1)]
%two different solutions - not surprising for a nonlinear equation...

%let's plot the graphs of s|-->u(1;s) the solution of iVp with u'(0)=s and u(0)=0
s=linspace(-4,40)';
xs1=zeros(size(s));
for k=1:100,
  X=lsode(@(x,t) [x(2);-x(1)*x(1)+f(t)],[0;s(k)],[0,1]);
  xs1(k)=X(2,1);
end  
plot(s,xs1,";u(1,s);",[s0,s1],[0,0],"*");pause(1);

%looks lie a parabol? let fit the square polynomial
p2=polyfit(s,xs1,2);
plot(s,xs1,s,polyval(p2,s),";poly;")
pause(1);
%so it isnot true - far from a parabole 
%perhaps a higher degree poly? take 4
p1=polyfit(s,xs1,4);
plot(s,xs1,";u(s);",s,polyval(p2,s),";square poly;",s,polyval(p1,s),";degree 4;")
