bagleytorvikequation(A,B,C)

function Y0 = bagleytorvikequation(A,B,C)
% 
% Sample function for solving Bagley-Torvik equation with zero initial conditions:
% A,B,C are coefficients of the Bagley-Torvik equation
% 
% Ay'' + By^(3/2) + Cy = F(t)
% 
% Reference: 
% I. Podlubny, Fractional Differential Equations, Academic Press, 
% San Diego, 1999, ISBN 0125588402. 
% (Correct typos: in Fig. 4 there on page 230 the results for A=B=C=1 are shown,
% replace the text "B=0.5" and "C=0.5" on page 231 with correct "B=1" and
% "C=1")

alpha=1.5; 

% Numerical solution:
h=0.075;
T=0:h:30;
N=30/h + 1;
M=zeros(N,N); % pre-allocate matrix for the system

% (1) Make the matrix for the alpha-th derivative:
Dalpha = ban(alpha,N,h);

% (2) Make the matrix for the second derivative:
D2 = ban(2,N,h);

% (3) Make the matrix for the entire equation:
M=A*D2 + B*Dalpha + C*eye(size(Dalpha));

% Make right-hand side:
F=8*(T<=1)';

% (3) Utilize zero initial conditions:
M = eliminator(N,[1 2])*M*eliminator(N, [1 2])';
F= eliminator(N,[1 2])*F;

% (4) Solve the system BY=F:
Y=M\F;

% (5) Pre-pend the zero values (those due to zero initial conditions)
Y0 = [0; 0; Y];

% Plot the solution:
figure(1)
plot(T,Y0,'g')
grid on