function [t,y]=bvp_df(A,q,Ba,Bb,r,a,b,l,h)
% Risolve il BVP lineare
% y(t)'=A(t,l)y(t)+q(t,l) a< t< b, con l parametro fissato
% Ba y(a)+ Bb y(b) = r
% con il metodo delle differenze finite basato sul metodo del punto di
% mezzo, con passo di discretizzazione h.
% In input: - le funzioni A(t,l) e q(t,l) che definisco l'eq differenziale
%           - le matrici Ba e Bb e il vettore r che definiscono le condizioni
%             al contorno
%           - gli estremi a e b dell'intervallo d'integrazione
%           - il passo h d'integrazione
% In output: - i vettori t e y, tali che y(:,n) e' l'approssimazione della
%              soluzione y in t(n) 
%
% Prova con Es. 7.2 Ascher Petzold: le function A.m e q.m definiscono A e
% q; la function bc.m restituisce in output B0, B1, e il vettore t
% Ad esempio provare:
% l=1; h=0.1; [B0, B1, b]=bc(l); [t,x]=bvp_df(@A, @q, B0, B1, b, 0, 1, l, h);
% Calcolare l'errore dell'approsimazione, sapendo che la soluzione esatta
% e'
% u=(exp(l*(t-1))+2*exp(2*l*(t-1))+exp(-l*t))/(2+exp(-l))+cos(pi*t);

t=[a:h:b];
m=length(r);
N=length(t)-1;
S=zeros(m,m,N);
R=zeros(m,m,N);
f=zeros(m*(N+1),1);
for n=1:N
    t1=t(n)+h/2;
    w=feval(A,t1,l);
    S(:,:,n)=-(eye(m)/h+0.5*w);
    R(:,:,n)=eye(m)/h-0.5*w;
    f((n-1)*m+1:n*m)=feval(q,t1,l);
end
f(N*m+1:(N+1)*m)=r;
H=spalloc(m*(N+1),m*(N+1),2*m^2*(N+1));
for n=1:N
    H((n-1)*m+1:n*m,(n-1)*m+1:n*m)=S(:,:,n);
    H((n-1)*m+1:n*m,n*m+1:(n+1)*m)=R(:,:,n);
end
H(N*m+1:(N+1)*m,1:m)=Ba;
H(N*m+1:(N+1)*m,N*m+1:(N+1)*m)=Bb;
w=H\f;
y=zeros(m,N+1);
for n=1:N+1
    y(:,n)=w((n-1)*m+1:n*m);
end

