function shoot_es3(l)
% Example 7.2 Ascher Petzold
% x'(t)=A(t)x(t)+q(t)
% Condizioni al contorno x1(0)=b1, x1(1)=b2; x2(1)=b3
% A(t) e q(t) dipendono dal parametro l
% Il metodo di shooting converge in una iterazione perche'
% il problema e' lineare
% 
% Prova l=1, 10, 20
% Shooting funziona con l=1;
% fallisce con l=10 (prova a diminuire AbsTol, RelTol)
% Differenze finite funziona molto bene (lo vedremo la prox settimana!)

A=zeros(3);
A(1,2)=1;A(2,3)=1;
A(3,:)=[-2*l^3, l^2, 2*l];

B0=zeros(3);
B0(1,1)=1;
B1=zeros(3);

B1=zeros(3); B1(2,1)=1; B1(3,2)=1;

% Le BC sono B0 x(0) + B1 x(1)= b
b1=(exp(-l)+exp(-l*2)+1)/(2+exp(-l))+1;
b2=0;
b3=(-l*exp(-l)+3*l)/(2+exp(-l));
b=[b1;b2;b3];

%Jacobiano
Q=B0+B1*expm(A);
cond(Q), eig(A)

options=[];

% Risolvo l'IVP con s=[0; 0; 0]
[t,x]=ode45(@es_bvp,[0,1],[0;0;0],options,l);

% B1*x(end,:)'-b e' il defect relativo a s=[0; 0; 0]
s=-inv(Q)*(B1*x(end,:)'-b);

% Newton converge in una iterazione, quindi s contiene
% i valori iniziali "esatti"
clear t

[t,y]=ode45(@es_bvp,[0,1],s,options,l);

% Soluzione esatta
u=(exp(l*(t-1))+exp(2*l*(t-1))+exp(-l*t))/(2+exp(-l))+cos(pi*t);

plot(t,y(:,1),'r',t,u,'g')

error=max(abs(y(:,1)-u))


function z=es_bvp(t,x,l)
A=zeros(3);
A(1,2)=1;
A(2,3)=1;
A(3,:)=[-2*l^3, l^2, 2*l];

q=zeros(3,1);
q(3)=sin(pi*t)*pi^3+ 2* l^3 *cos(pi*t) + l^2*sin(pi*t)*pi+2*l*cos(pi*t)* pi^2;

z=A*x+q;

