% Script for solving Bagley-Torvik equation with zero initial conditions:

clear all

alpha=1.5;

% Numerical solution:
h=0.075;
T=0:h:30;
N=30/h + 1;
B=zeros(N,N);

% (1) Make the main matrix B1_{N}^{alpha}:
bc=fliplr(bcrecur(alpha,N-1));
for k=1:N
   B1(k,1:k)=bc((N-k+1):N);
end

% (2) Make the main matrix B2_{N}^{2}:
alpha=2;
bc=fliplr(bcrecur(alpha,N-1));
for k=1:N
   B2(k,1:k)=bc((N-k+1):N);
end

% (3) Make the matrix for the entire equation:
B=h^(-2)*B2 + h^(-1.5)*B1 + eye(size(B1));

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

% (3) Use zero initial conditions:
F(1)=0;
F(2)=0;
B(1,1)=1;
B(2,1)=-1/h; B(2,2)=1/h;

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

% Plot the solution:
figure(1)
plot(T,Y,'r')



