//
// funzioni di servizio
//

function [mesh,h] = mesh1D (a,b,N)
//costruzione della mesh
h = (b-a)/N;
mesh = zeros(1,N+1);
for i=1:N+1
    mesh(i) = a+(i-1)*h;
end
endfunction


function [Ki,Hi,Mi] = matrices1D (h)
// Definizione delle matrici locali
Ki = [1 -1; -1 1]/h;
Hi= [-1 1; -1 1]/2;
Mi = h*[2 1; 1 2]/6;
endfunction


function [K,H,M] = assemble1D (Ki,Hi,Mi,N)
// Definizione delle matrici globali

 K = zeros(N+1,N+1);
 H = zeros(N+1,N+1);
 M = zeros(N+1,N+1);

for i=1:N
    K(i:i+1,i:i+1) = K(i:i+1,i:i+1) + Ki;
    H(i:i+1,i:i+1) = H(i:i+1,i:i+1) + Hi;
    M(i:i+1,i:i+1) = M(i:i+1,i:i+1) + Mi;
end
endfunction


function [fi] = localf (mesh,h,beta,gamma,u_a,u_b,a,b)
//
// vettore con i valori di f 
//(opportunamente corretta se non omogeneo) sulla mesh
//
  retta=u_a+(mesh-a)*(u_b-u_a)/(b-a);
  der_retta=mesh*0+(u_b-u_a)/(b-a);
  fi = (h/2)*(f(mesh)-beta*der_retta-gamma*retta)';
endfunction


function [rhs] = assemblef (fi,N)
// Costruzione del termine noto f

rhs=zeros(N+1,1);
for i=1:N
    rhs(i:i+1,1) = rhs(i:i+1,1) + fi(i:i+1,1);
end
endfunction

//---------------------------------------------------------
function [value] = uhP1(t, uh, n,a,b)
// implement the piecewise linear function u_h in the point t
    h=(b-a)/n;
for i=1:length(t)
    k=floor((t(i)-a)/h)+1;
    if abs(t(i)-b)<.000001
        k=n;
    end
    uelem(1) = uh(k);
    uelem(2) = uh(k+1);
    t_loc=(t(i)-(a+(k-1)*h))/h;
    [b1, b2]= basisP1(t_loc);
    value(i) = (uelem)' * ([b1, b2])';
end
endfunction
//----------------------------------------------------------
function [b1, b2, bx1, bx2]= basisP1(t)
//Evaluation of P1 basis defined in [0,1]
b1 = (1-t);
b2 = t;
// Evaluation of the derivatives
bx1=-1;
bx2=1;
endfunction
//--------------------------------------------------------------
function [err] = Errore_L2 (a,b,h,uh,N)
// calcolo dell'errore in norma L2  (con trapezi) 
h_tab=h/2
tab_err=[a:h_tab:b];
uesa=usol(tab_err);
ucalc=uhP1(tab_err,uh,N,a,b);
err=0;
for i=1:max(size(tab_err))-1
    err= err+((uesa(i)-ucalc(i))^2+(uesa(i+1)-ucalc(i+1))^2)*h_tab/2;
end
err=err^(.5);
endfunction
//-----------------------------------------------------------------
// FUNZIONI da CAMBIARE a SECONDA del PROBLEMA
//--------------------------------------------------------------------------------------
function [value] = f(t)
// valutazione del termine noto f
  value=-4*%pi*cos(2*%pi*t)+(4*(%pi^2)+1)*t.*sin(2*%pi*t);
endfunction
//----------------------------------------------------------------
function [value] = usol(t)
// valutazione della corrispondente soluzione esatta
// alfa=1, beta=0, gamma=1, 
 value = t.*sin(2*%pi*t);
endfunction 

//---------------------------------------------------------------    
//Soluzione problema ellittico omogeneo unidimensionale
//----------------------------------------------------------------
// -alfa u''+beta u' +gamma u=f in (a,b)
// u(a)=u_a(b)=u_b;
//----------------------------------------------------------------
//
// si immagina di conoscere la soluzione esatta
//
//-----------------------------------------------------------------------
function [uh,h,mesh,errore,retta] = elliptic1D(a,b,u_a,u_b,N,alfa,beta,gamma)
//
// costruzione della mesh
[mesh,h] = mesh1D(a,b,N);
//
// costruzione della matrice del sistema linerae 
[Ki,Hi,Mi] = matrices1D(h);
[K,H,M] = assemble1D(Ki,Hi,Mi,N);
Atilde = alfa*K +beta*H+gamma*M;
fine=max(size(Atilde));
A = Atilde(2:fine-1,2:fine-1);
num_cond=cond(A);
//
// costruzione del termine noto f
fi = localf(mesh,h,beta,gamma,u_a,u_b,a,b);
rhs = assemblef (fi,N);
rhs = rhs(2:fine-1,1);
//
//  risoluzione del sistema lineare
uh = A\rhs;
// imposizione delle condizioni al bordo  Dirichlet 
uh = [0; uh; 0];
retta=(u_a+(mesh-a)*(u_b-u_a)/(b-a))';
uh=uh+retta;
//
// Grafico della soluzione calcolata e di quella esatta
tt=[a:h/10:b];
plot(tt,uhP1(tt,uh,N,a,b),tt,usol(tt), 'r')
//
// calcolo dell'errore in norma L2 (con trapezi) se la soluzione esatta e' nota
errore=Errore_L2(a,b,h,uh,N);
//
endfunction
