

disp("linear shooting");
graphics_toolkit("gnuplot");
%solving -u''+u=0 u(0)=u(T)=1

vf=@(x,t) [x(2);x(1)]

T=1;
for k=1:9,
  s=0;
  X=lsode(vf,[1;s],[0,T]);
  A=X(2,1);
  s=1;
  X=lsode(vf,[1;s],[0,T]);
  B=X(2,1)-A;
  s=(1-A)/B
  t=linspace(0,T,200);
  X=lsode(vf,[1;s],t);
  [T,X(end,1)-1]
  plot(t,X(:,1),[";Shooting u''=u u(0)=u(",num2str(T),")=1;"],[0;T],[1,1],"*");pause(2);
  T+=2;
end

%nonlinear shooting
%u''=3/2 u^2 u(0)=4 u(1)=1
global vf
vf=@(x,t) [x(2);3*x(1)*x(1)*0.5];

function y=G(s)
  global vf 
  X=lsode(vf,[4;s],[0,1]);
  y=X(2,1)-1;
end  

t=linspace(0,1,300);
s1=fsolve(@G,-1)
X1=lsode(vf,[4;s1],t);
s2=fsolve(@G,-30)
X2=lsode(vf,[4;s2],t);
plot(t,X1(:,1),t,X2(:,1),";Shooting u''=(3/2)u^2 u(0)=4,u(1)=1;");pause(2);


%FDM for a linear BVP -u''+cu=f in (a,b) u(a)=alpha u(b)=beta

function [u,x]=fdmbvpsolve(f,a,b,N,c=0,al=0,be=0)
%a function solving BVP   -u''+cu=f in (a,b) u(a)=alpha u(b)=beta
%by FDM on the uniform grid of N+! points in [a,b] i.e.
%we want to find [u(0),..,u(N)] s.t. u(k) \approx u(x(k)) for x(k)=a+k*h h=(b-a)/N
%INPUT
%f - a function handle to the rhs function f
%a,b- the interval endswitch
%N - no of gridpoints minus one
%c - coefficient of ODE (default =0)
%al be - boundary values at the ends (default =0)
%OUTPUT
%u - the solution vectorize
%x - vector of gridpoints
%EXAMPLE
%
%a=0;b=1;c=1;sol=@sin; f=@(x) (1+c)*sin(x);al=sol(a);be=sol(b);N=20;[u,x]=fdmbvpsolve(f,a,b,N,c,al,be);
%plot(x,u,x,sol(x));pause(1);plot(x,u-sol(x),";error;");

x=linspace(a,b,N+1)';
h=(b-a)/N;
C=A=sparse(N+1,N+1);
A(2:N,2:N)=sparse(diag(2*ones(N-1,1))) - sparse(diag(ones(N-2,1),-1))-sparse(diag(ones(N-2,1),1));
A(2,1)=A(N,N+1)=-1;
A*=(1/(h*h));
A(1,1)=A(N+1,N+1)=1.0;
C(2:N,2:N)=sparse(c*eye(N-1,N-1));

F=[al;f(x(2:N));be];
u=(A+C)\F;
end

%solving u''=u u(0)=u(T)=1

T=1; f=@(x) zeros(size(x));
for k=1:9,
 N=20*T;
 [u,x]=fdmbvpsolve(f,0,T,N,c=1,1,1);
 plot(x,u,[";FDM u''=u u(0)=u(",num2str(T),")=1;"]);pause(2);
 T+=3;
end
