function rt = restrid(n,alp)
% n   = number of intervals
% alp = const > 0

a=zeros(n+1);
b=zeros(n+1);
c=zeros(n+1);
r=zeros(n+1);
    
if (alp<=0) 
    fprintf('Alpha must be > 0 !!! \n');
    return;
end;

a(2:n) = -1.0;
b(2:n) =  2.0;
c(2:n) = -1.0;
for i=2:n
    r(i) = - f_right( (i-1)/n , alp ) / n^2 ;
end;

fprintf('\n We suppose zero values at boundary \n\n');

c(1)=0.0;
b(1)=1.0;   b(n+1)=1.0;
r(1)=0.0;   r(n+1) = 0.0;

u = tridiag(n+1, a,b,c,r);

maxr=0.0;
maxh=0.0;

chk  = zeros(size(u));
xval = zeros(size(u));
for i=1:n+1
    xval(i) = (i-1)/n;
    chk(i) = (1.0+xval(i)*(f_right(1.0,alp)-1.0)-f_right(xval(i),alp)) / alp^2;
    fprintf('%5d %9.4f:   y(i)= %+12.7e,  check= %+12.7e\n',i,xval(i),u(i),chk(i));
    maxr = max( maxr, abs(u(i)-chk(i)) );
    maxh = max( maxh, abs(u(i))     ) ;
end

fprintf('\n Maximum difference / maximum of result = %12.7e \n\n',maxr/maxh);
    
figure;
plot(xval,chk,'k');
hold on;
plot(xval,u,'ro')
hold off

return;


%%%%%%%%%%%%

function fr = f_right(x,alp);
    fr = cos(alp*x);

 
%%%%%%%%%%%%

function u = tridiag(n, a,b,c,r);
% Find solution u of a tridiagonal system of equations
%
% Vector a is below the diagonal, a(1) unused
% Vector b is on the diagonal
% Vector c is above the diagonal, c(n) unused 
% Vector r is the right-hand side    

u     = zeros(n);
gamma = zeros(n);

if (b(1)==0.0)     
    fprintf('TRIDIAG: b(1)=0 \n');
    pause
end;
beta  = b(1);
u(1) = r(1) / beta;
for j=2:n
    gamma(j) = c(j-1)/beta;
    beta = b(j) - a(j)*gamma(j);
    if (beta==0.0)
        fprintf('TRIDIAG: beta=0 \n');
        pause
    end;
    u(j) = (r(j)-a(j)*u(j-1)) / beta;
end;
for j=n-1:-1:1
    u(j) = u(j) - gamma(j+1)*u(j+1);
end;


%%%%%%%%%%%%

